COMPUTATIONAL  STUDY  OF  LEUKOCYTE  RHEOLOGY 
BASED  ON  A  MULTILAYER  FLUID  MODEL 


■  By 
HENG-CHUAN  KAN 


A  DISSERTATION  PRESENTED  TO  THE  GRADUATE  SCHOOL 
OF  THE  UNIVERSITY  OF  FLORIDA  IN  PARTIAL  FULFILLMENT 
OF  THE  REQUIREMENTS  FOR  THE  DEGREE  OF 
DOCTOR  OF  PHILOSOPHY 

UNIVERSITY  OF  FLORIDA 
1997 


ACKNOWLEDGMENTS 

I  would  like  to  express  my  sincere  thanks  toward  my  advisor.  Dr.  Roger  Tran-Son-Tay, 
for  his  expert  guidance,  encouragement,  and  financial  support  during  the  course  of  my 
graduate  studies.  He  has  been  a  constant  source  of  inspiration  for  me  both  professionally  and 
philosophically.  The  same  gratitude  is  also  given  to  my  cochair.  Dr.  Wei  Shyy,  for  his 
understanding  and  help  in  so  many  ways. 

Very  special  thanks  go  to  Dr.  H.  S.  Udaykumar  for  his  tremendous  help  and  enthusiasm 
in  this  work.  In  spite  of  spending  the  countless  and  frustrated  hours  in  front  of  a  variety  of 
computers  with  me  during  the  early  debugging  process,  his  patience  and  belief  kept  this  work 
going. 

I  would  like  to  thank  the  members  of  my  supervisory  committee.  Dr.  Renwei  Mei,  Dr. 
Richard  B.  Dickinson  and  Dr.  Nicolae  D.  Cristescu,  for  their  fruitful  advices. 

I  wish  to  thank  my  colleagues  at  the  Laboratory  of  Cellular  Mechanics  and  Biorheology 
for  their  direct  and  indirect  contributions  and  entertaining  discussions.  In  particular,  I  thank 
Emmanuelle,  David,  Neal,  Jie,  Philou,  and  Philippe  for  their  friendship  and  assistance. 

Words  are  not  enough  to  describe  my  sincere  gratitude  to  my  parents,  Pao-Hsiang  Wang 
and  Che-Chang  Kan,  for  their  persistent  support  and  love.  Their  emphasis  on  higher 
education  for  their  children  pave  the  opportunity  for  me  to  complete  this  work.  I  also  wish  to 
thank  my  brothers,  Po-Wen  Kan  and  Feng-Jung  Kan,  for  their  continuous  encouragement 
throughout  these  years. 

Finally,  I  reserve  my  deepest  appreciation  to  my  wife,  Min-Hui.  I  will  be  forever 
indebted  to  her  endless  sacrifices  and  love.  To  her,  I  dedicate  this  dissertation. 


ii 


TABLE  OF  CONTENTS 


page 

ACKNOWLEDGMENTS    ii 

ABSTRACT    v 

L     INTRODUCTION  AND  BACKGROUND   1 

LL  Motivation  and  Objectives   1 

1.2.  Importance  of  Compound  Drops  Dynamics    2 

1.3.  Significance  of  Leukocyte  Rheology   2 

1 .4.  Observed  Rheological  Behavior  of  Leukocytes    4 

1.5.  Current  Status  of  Leukocyte  Modeling   6 

1.6.  Summary    9 

2  THEORETICAL  ANALYSIS   12 

2.1.  Background    12 

2.2.  ELAFINT:  A  Mixed  Eulerian-Lagrangian  Solution  Algorithm   15 

2.2.1.  Characteristics  of  the  Numerical  Methodology   15 

2.2.2.  Interface  Tracking  Procedure   16 

2.2.3.  Formulation  of  the  Pressure-Based  Flow  Solver    18 

2.2.4.  Discretization  Procedure  of  the  Flow  Field  Equations    20 

2.3.  Momentum  Interpolation   23 

2.4.  Interactions  between  Interfacial  Physics  and  Flow  Field  Variables    28 

2.4. 1 .  Cut-Cell  Procedures  for  Solid-Liquid  Boundaries   28 

2.4.2.  Immersed  Boundary  Technique  for  Fluid-Fluid  Interfaces   31 

2.5.  Mass  Conservation  in  Individual  Phases   36 

2.6.  Summary    39 

3  ASSESSMENT  OF  COMPUTATIONAL  TECHNIQUE   40 

3.1.  Flows  With  Fixed  Arbitrary-Shaped  Internal  Boundaries    40 

3.2.  Moving  Boundary  Flow  Problems   47 

3.2. 1 .  Deformation  and  Recovery  of  Viscous  Drops  in  Extensional  Flows  47 

iii 


3.2.2.  Deformation  of  Drops  in  Extensional  Flows  with  Inertia  Effects  .  58 

3.2.3.  Deformation  of  Drops  in  Constricted  Tubes   65 

3.3.  Summary    71 

4  DYNAMICS  OF  COMPOUND  DROPS    72 

4.1.  Problem  Statement  and  Mathematical  Formulation    72 

4. 1 . 1 .  Computational  Setup   76 

4. 1 .2.  Selection  of  Model  Parameters    77 

4.2.  Deformation  Analysis    77 

4.2.1.  Steady-State  Configuration  (  Ca  <  Cacr )    77 

4.2.2.  Unsteady  Response  (  Ca  >  Cacr )   79 

4.3.  Recovery  Analysis   81 

4.3.1.  Effect  of  the  Capillary  Number   81 

4.3.2.  Effect  of  the  Core  Viscosity   84 

4.3.3.  Comparison  of  Simple  and  Compound  Drops   86 

4.4.  Summary    88 

5  LEUKOCYTE  MODELING   89 

5.1.  Choice  of  Model  Parameters   89 

5.2.  Effect  of  a  Viscous  Cytoplasmic  Layer   90 

5.3.  Effects  of  Nucleus  Deformation  and  Time  Scales   93 

5.3.1.  Incompatible  Nucleus  and  Cytoplasm  Time  Scales   93 

5.3.2.  Compatible  Nucleus  and  Cytoplasm  Time  Scales   96 

5.4.  Discussion  And  Implications  to  Leukocyte  Rheology    102 

6  SUMMARY  AND  FUTURE  RESEARCH   105 

6.1.  Conclusion   105 

6.2.  Future  Research   106 

6.2.1.  Experimental  Work    107 

6.2.2.  Numerical  Algorithm   107 

REFERENCES    109 

BIOGRAPHICAL  SKETCH   117 


iv 


Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirement  for  the  Degree  of  Doctor  of  Philosophy 

COMPUTATIONAL  STUDY  OF  LEUKOCYTE  RHEOLOGY 
BASED  ON  A  MULTILAYER  FLUID  MODEL 

By 

Heng-Chuan  Kan 
May  1997 

Chairperson:  Dr.  Roger  Tran-Son-Tay 

Major  Department:  Aerospace  Engineering,  Mechanics  and  Engineering  Science 

Knowledge  of  the  rheological  properties  of  leukocytes  (white  blood  cells)  is  essential 
not  only  for  the  comprehension  of  microcirculatory  flow  dynamics,  but  also  for  the  under- 
standing of  their  functions  and  behavior  in  health  and  disease.  These  properties,  together 
with  the  leukocyte's  structural  characteristics,  determine  the  deformability  of  the  cell,  espe- 
cially during  large  deformations,  such  as  those  involved  in  their  release  from  the  bone  mar- 
row or  extravasation  into  the  interstitium.  Existing  viscoelastic  solid,  Newtonian  and  non- 
Newtonian  liquid  drop  models  used  to  calculate  the  cell  rheological  properties  do  not 
adequately  explain  reported  experimental  observation.  This  is  because  such  models  do  not 
take  into  account  the  morphology  of  the  cell  and  its  multi-layer  structure. 

In  this  work,  the  deformation  and  recovery  of  a  compound  liquid  drop  (three-layer  fluid 
model)  are  investigated  by  means  of  a  mixed  Eulerian-Lagrangian  computational  methodol- 
ogy which  allows  large  viscosity  and  capillarity  differences  between  layers.  Analysis  of  the 
deformation  process  confirms  published  results  and  clearly  indicates  that  the  history  of  the 
deformation  is  stored  in  the  compound  drop  in  terms  of  the  shape  of  its  components.  The 
recovery  process  is  found  to  provide  information  on  the  hydrodynamics  of  compound  drops 


that  is  not  supplied  by  the  deformation  analysis.  It  is  discovered  that  a  compound  drop  does 
not  recover  like  a  single  phase  Newtonian  liquid  drop  except  when  the  core  is  sufficiently 
deformed  and  its  material  properties  are  such  that  the  core  deformation/recovery  time  scale 
is  compatible  with  that  of  the  shell  layer.  A  difference  in  the  core  and  shell  layer  time  scales 
causes  an  initial  rapid  recoil  of  the  drop  during  which  the  shell  fluid  is  the  sole  participant 
in  the  hydrodynamics,  followed  by  a  slower  relaxation  period  during  which  the  core  and 
shell  layer  couple.  The  large  deformation  and  subsequent  recovery  processes  of  compound 
drops  also  provide  new  insights  into  the  complicated  behavior  of  leukocytes.  This  three-lay- 
er fluid  model  describes  the  morphology  of  leukocytes  better  than  existing  ones.  It  consists 
of  an  outer  interface  (cell  membrane),  containing  a  shell  layer  (cytoplasm),  and  a  core 
(nucleus). 

Based  on  the  compound  drop  dynamics,  it  is  concluded  that  (1)  the  nucleus  plays  an 
essential  role  in  the  cell  response,  and  (2)  the  energy  stored  in  the  deformed  interfaces  can 
be  important.  The  present  investigation  indicates  that  unless  the  nucleus  and  its  deformation 
are  included  in  the  analyses,  neither  Newtonian  nor  non-Newtonian  drop  models  for  leuko- 
cytes can  yield  a  reliable  picture  of  the  hydrodynamics  of  such  cells. 


vi 


CHAPTER  1 
INTRODUCTION  AND  BACKGROUND 

1.1.  Motivation  and  Objectives 

Knowledge  of  the  rheological  properties  of  leukocytes  (white  blood  cells)  is  essential 
not  only  for  the  comprehension  of  microcirculatory  flow  dynamics,  but  also  for  the  under- 
standing of  their  functions  and  behavior  in  health  and  disease.  These  properties,  together 
with  the  leukocyte's  structural  characteristics,  determine  the  deformability  of  the  cell,  espe- 
cially during  large  deformations,  such  as  those  involved  in  their  release  from  the  bone  mar- 
row or  extravasation  into  the  interstitium.  This  deformability  is  also  essential  in  the  leuko- 
cyte's response  to  disease/infection  in  humans  because  it  determines  the  ability  of  the  cell 
to  flow  and  deform  in  capillaries  and/or  migrate  in  tissue.  Therefore,  obtaining  the  correct 
mechanical  description  for  the  leukocytes  has  been  a  focus  of  research  in  the  recent  years. 

Existing  viscoelastic  solid,  Newtonian  and  non-Newtonian  liquid  drop  models  used  to 
calculate  the  cell  rheological  properties  do  not  adequately  explain  reported  experimental  ob- 
servation. This  is  because  such  models  do  not  take  into  account  the  morphology  of  the  cell 
and  its  multilayer  structure. 

The  three-layer  (or  compound  drop)  model  has  been  suggested  to  be  a  good  model  for 
the  leukocytes  (Skalak  et  al.,  1990;  Dong  et  al.,  1991;  Hochmuth  et  al.,  1993).  This  model 
gives  a  better  representation  of  the  morphology  of  the  leukocyte  and  its  interior  over  the  ex- 
isting ones.  It  consists  of  a  cell  membrane  (outer  interface),  containing  a  cytoplasmic  layer 
(shell),  and  a  nucleus  (core  and  inner  interface).  However,  the  three-layer  model  has  not 
hitherto  been  approached  analytically  due  to  the  limitation  of  available  tools. 

The  objectives  of  this  work  are:  (1)  to  develop  a  general  computational  methodology 
for  investigating  the  rheological  behavior  of  multilayer  drop  systems,  (2)  to  test  the  validity 
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and  accuracy  of  the  algorithm,  (3)  to  study  the  hydrodynamics  of  compound  drops,  and  (4) 
to  provide  a  better  understanding  of  leukocyte  rheology. 

1.2.  Importance  of  Compound  Drops  Dynamics 

The  deformation  and  subsequent  relaxation  of  a  compound  liquid  drop  is  a  fundamental 
fluid  mechanics  problem  of  medical  importance.  The  underlying  idea  that  makes  the  ap- 
plication of  compound  drops  attractive  and  has  spawned  several  new  technologies  (Jackson 
and  Lee,  1991;  Powell  and  Mahalingam,  1992;  Batich  and  Vaghefi,  1994;  Goeletal.,  1994; 
Martin  and  Parthasarathy,  1995;  Hassan  et  al.,  1996;  Thies,  1996)  is  that  a  core  of  active  in- 
gredient (whether  solid,  liquid  or  gas)  can  be  protected  by  enclosing  it  within  a  shell  of  a 
second  material.  The  active  ingredient,  which  must  not  of  course  attack  the  chosen  encapsu- 
lating material,  is  released  at  a  determined  rate,  time  and  location  to  achieve  certain  localized 
effects.  These  include  controlled  release  of  drugs  (Tai  and  Sun,  1993),  flavors  (Jackson  and 
Lee,  1991;  Gorski,  1994),  enzymes  (Hassan  et  al.,  1996)  or  latent  heat  (Goel  et  al.,  1994; 
Mulligan  et  al.,  1996). 

While  many  aspects  of  compound  drops  have  been  studied  in  the  literature  (Johnson 
and  Sadhal,  1985),  including  in  particular  deformation  in  an  extensional  flow  (Stone  and 
Leal,  1990)  and  stability  (Landman,  1985;  Oguz  and  Sadhal,  1987),  the  range  of  conditions 
under  which  compound  drop  dynamics  has  been  investigated  does  not  span  the  parameters 
of  interest  in  this  work. 

1.3.  Significance  of  Leukocyte  Rheology 

The  persistence  and  widespread  occurrence  of  diseases  associated  with  defects  in  the 
circulation  of  blood  has  largely  provided  the  impetus  for  the  study  of  blood  rheology.  Human 
blood  is  a  suspension  of  erythrocytes  (red  blood  cells),  leukocytes  (white  blood  cells)  and 
platelets  in  plasma.  Generally  speaking,  blood  is  considered  as  a  Newtonian  fluid  in  large 
vessels  at  shear  rates  greater  than  100  s~' .  When  the  shear  rate  is  less  than  10  s~' ,  blood  be- 
haves like  a  Casson  fluid  (Cokelet  et  al.,  1963).  However,  for  blood  vessels  less  than  500 
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\im  in  diameter,  blood  cannot  be  considered  as  a  continuous  fluid.  Blood  must  be  treated 
as  a  two  phase  fluid  (core  and  plasma  layers)  in  vessels  with  diameters  between  30  to  500 
^m.  Several  effects  (e.  g.,  Fahraeus,  and  Fahraeus-Lindqvist)  occur  in  that  diameter  range 
and  become  stronger  as  the  vessel  diameter  decreases  (Fung,  1993).  For  blood  vessels  less 
than  20  ~  30  \im  in  diameter,  the  presence  of  individual  cells  cannot  be  neglected.  The  flow 
is  creeping  and  the  rheological  properties  of  individual  blood  cells  have  to  be  taken  into  ac- 
count. 

The  rheological  behavior  of  erythrocytes  has  been  intensively  studied  and  found  to  be 
adequately  modeled  by  a  Newtonian  fluid  enclosed  in  a  viscoelastic  membrane  (Evans  and 
Hochmuth,  1976;  Evans  and  Skalak,  1980).  On  the  other  hand,  the  search  for  a  suitable  mod- 
el to  describe  the  rheological  behavior  of  leukocytes  continues. 

Rheological  studies  of  leukocytes  have  been  neglected  because  they  occupy  less  than 
1%  of  the  total  volume  content  of  blood.  However,  since  leukocytes  have  found  to  have  a 
larger  volume  and  a  lower  deformability  than  red  blood  cells,  the  number  of  studies  has 
steadily  increased.  Leukocytes,  because  of  their  relative  size  and  rigidity,  play  a  significant 
role  on  the  overall  circulatory  flow.  They  can  affect  in  particular  the  apparent  viscosity  of 
blood  and  therefore  the  pressure  drops  encountered  across  vessels  and  capillaries  (Pries  et 
al.,  1994).  In  addition,  the  ability  of  the  leukocytes  to  flow  and  deform  in  capillaries  and/or 
migrate  in  tissue  is  essential  in  the  response  to  disease/infection  in  humans.  Therefore, 
knowledge  of  the  rheological  properties  of  the  leukocytes  is  critical  for  understanding  micro- 
circulatory  flow  dynamics  and  its  effects  in  health  and  disease. 

Given  the  microscopic  size  of  the  blood  cells  (about  8—10  [im  in  diameter),  experimen- 
tal efforts  face  obvious  hurdles.  However,  micro-manipulation  techniques  have  been  suc- 
cessfully developed  not  only  for  deforming  but  also  for  evaluating  mechanical  and  chemical 
effects  on  blood  cells  (Schmid-Schonbein  et  al.,  1981;  Evans  and  Kukan,  1984;  Needham 
and  Hochmuth,  1990;  Usami  et  al.,  1992;  Skierczynski  et  al.,  1993;  Tran-Son-Tay  et  al., 
1994a).  One  common  way  of  determining  the  rheological  properties  of  the  leukocyte  has 
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been  to  deform  the  cell  by  aspiration  into  a  micropipet  (3  ~  6  |xm  in  diameter)  and  then  to 
eject  the  cell  and  observe  the  recovery  characteristics.  The  recovery  time  and  path  then  pro- 
vide an  estimate  of  the  "apparent  viscosity"  of  the  cell.  The  apparent  viscosity  is  an  average 
rheological  property  of  the  whole  cell,  including  all  the  inclusions. 

1 .4.  Observed  Rheological  Behavior  of  Leukocytes 

The  most  populous  type  of  leukocytes  in  human  blood  is  the  neutrophil.  The  neutrophil 
contains  a  relatively  small  segmented  nucleus,  occupying  about  21%  of  its  volume,  and  a 
cytoplasm  that  is  rich  in  granules.  In  contrast,  in  a  lymphocyte,  another  type  of  leukocyte, 
the  nucleus  occupies  approximately  44%  of  the  cells  (Schmid-Schonbein  et  al.,  1980).  A 
neutrophil  is  spherical  in  a  passive,  non-reactive  state.  The  membrane  of  the  neutrophil  is 
highly  ruffled,  giving  it  an  excess  surface  area  (estimated  at  110%—  120%  of  the  unde- 
formed  sphere),  and  permits  deformation  at  constant  volume. 

In  spite  of  the  intensive  research  effort,  to  date,  a  variety  of  behaviors  of  the  leukocyte 
cannot  be  explained  clearly.  For  example,  when  a  neutrophil  is  aspirated  into  a  micropipet 
it  flows  into  the  pipet  when  the  aspiration  pressure  exceeds  a  value  called  the  critical  suction 
pressure.  Neutrophils  rapidly  aspirated  into  a  micropipet  with  a  diameter  smaller  than  the 
cell  diameter  enter  at  an  initial  velocity  much  faster  than  the  steady-state  velocity.  Also,  as 
the  rate  of  deformation  increases,  the  cell  viscosity  is  found  to  decrease,  which  has  been  spec- 
ulated to  be  due  to  a  shear-thinning  effect  (Needham  and  Hochmuth,  1990;  Tsai  et  al.,  1993). 
Needham  and  Hochmuth  (1990)  observe  that  after  the  cell  enters  into  the  pipet,  the  subse- 
quent deformation  occurs  at  a  constant  viscosity,  independent  of  the  aspiration  pressure. 
However,  as  the  excess  aspiration  pressure  increases,  the  apparent  viscosity  of  the  cell  de- 
creases. Furthermore,  following  small  but  rapid  deformations  (in  less  than  0.3  seconds)  of 
the  cell,  a  rapid  initial  recovery  (around  0.3  seconds)  is  obtained,  followed  by  a  slower  phase 
(Dong  et  al.,  1991).  On  the  other  hand,  when  the  neutrophil  is  aspirated  and  held  inside  the 
pipet  for  a  short  time  period  of  5  seconds  or  less,  and  then  expelled,  the  cell  recovers  faster 
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Authors 

Experiments 

Cell  Apparent  Viscosity 
Operating  conditions  ^(Poise) 

Dong  etal.  (1988) 

Aspiration 
ixeL-uvciy 

Small  deformation 

300 

Evans  and  Yeung 
(1989) 

Aspiration 

Large  deformation 

1000-2000 

Needham  and 
Hochmuth  (1990) 

Aspiration 

High  deformation  rate 
Low  deformation  rate 

1000 
2000 

Dong  etal.  (1991) 

Aspiration 
Recovery 

High  deformation  rate 
Large  deformation 

O  AA    O  AAA 

zlHJ— zUUU 

Tran-Son-Tay 
etal.  (1991) 

Dong  and  Skalak 
(1992) 

Recovery 
Aspiration 

Large  deformation 
(short  holding  time) 

Tyarpe  deformation 
(long  holding  time) 

Large  deformation 

800-1500 
1000-2000 
50-600 

Hochmuth  et  al. 
(1993) 

Tsai  et  al.  (1993) 

Recovery 
Aspiration 

Small  deformation 

Low  aspiration  pressure 
High  aspiration  pressure 

600 

5000 
500 

Waugh  and  Tsai 
(1994) 

Aspiration 

High  deformation  rate 

1000-1500 

Table  I.  A  list  of  the  published  values  of  neutrophil  cytoplasmic  viscosity 
based  on  a  fluid  model  under  different  experimental  conditions. 


than  when  held  for  a  longer  time  with  an  initial  rapid  recoil  (Tran-Son-Tay  et  al.,  1991). 
Finally,  neutrophils  that  are  only  slightly  deformed  in  a  pipet  before  ejection,  yield  much 
smaller  values  of  apparent  viscosity  (Hochmuth  et  al.,  1993)  than  cells  that  are  significantly 
deformed.  This  is  analog  to  the  finding  that  the  cell  appears  to  be  more  viscous  when  it  flows 
into  smaller  diameter  pipets  (Evans  and  Yeung,  1989). 

It  appears  that  the  rate  of  deformation,  the  extent  of  deformation,  the  excess  aspiration 
pressure  and  the  holding  time  all  affect  the  cell  recovery  characteristics,  resulting  in  different 
apparent  viscosities.  Table  I  summarizes  the  apparent  viscosity  values  for  neutrophils  re- 
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ported  by  various  groups.  These  results  obtained  from  the  aspiration  (flow  into  a  pipet)  and 
recovery  experiments.  Therefore,  these  curious  properties  of  the  cell  have  yet  to  be  fully 
explained  by  existing  theoretical  models  based  on  the  Newtonian  or  the  non-Newtonian 
framework. 

1 .5.  Current  Status  of  Leukocyte  Modeling 

Much  effort  has  been  done  over  the  years  on  obtaining  the  correct  rheological  model 
for  leukocytes.  In  earlier  studies,  Bagge  et  al.  (1977)  show  that  granulocytes  deform  as  they 
flow  down  a  tapered  pipet.  Transit  times  down  the  tube  were  on  the  order  of  seconds,  for 
pressure  drops  of  magnitudes  found  in-vivo.  When  expelled  from  the  pipet,  the  granulo- 
cytes recover  their  original  spherical  shape  in  times  ranging  from  20  to  80  seconds.  In  the 
work  of  Schmid-Schonbein  et  al.  (1981),  the  leukocyte  is  modeled  as  a  simple  viscoelastic 
solid  with  a  Maxwell  element  in  parallel  with  an  elastic  element,  as  shown  in  Figure  1.1. 


Figure  1.1.  Illustration  of  a  standard  viscoelastic  solid  model  for  the  leukocyte. 
The  model  constants  \i  and  k's  represent  the  viscous  and  elastic  coefficients. 

The  elastic  element  is  thought  to  impart  shape  memory,  which  enables  the  cell  to  recoil  to 
the  resting  shape.  However,  Evans  and  Kukan  (1984)  show  that  the  deformation  of  a  leuko- 
cyte into  a  pipet  is  a  continuous  flow  process  with  no  approach  to  a  static  deformation  limit. 
It  is  therefore  suggested  that  the  leukocyte  be  modeled  as  a  Newtonian  liquid  drop  with  a 
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Figure  1.2.  Illustration  of  a  two-layer  Newtonian  liquid  model  with 
constant  cortical  tension  for  the  leukocyte,  y  represents  the  isotropic  sur- 
face tension  coefficient. 

prestressed  cortical  tension  (Yeung  and  Evans,  1989),  as  shown  in  Figure  1.2,  rather  than  a 
viscoelastic  solid.  The  shape  memory  then  would  be  induced  in  much  the  same  way  as  in 
drops,  by  ascribing  a  cortical  tension  (corresponding  to  surface  tension  in  drops)  to  the  mem- 
brane enclosing  the  cytoplasm.  Based  on  this  representation,  Yeung  and  Evans  (1989)  simu- 
late the  flow  of  a  neutrophil  into  a  pipet  by  imposing  a  constant  "ring"  force  at  the  pipet  tip 
and  a  constant  orifice  pressure  inside  the  pipet  to  satisfy  the  free  slip  condition. 

Although  the  liquid-like  behavior  of  neutrophils  is  widely  accepted,  the  precise  nature 
of  this  hquid  is  not  established.  In  their  experiments,  Needham  and  Hochmuth  (1990)  ob- 
served that  the  cell  exhibits  a  non-Newtonian  shear-thinning  behavior  when  flowing  into 
a  pipet  at  high  rates  of  deformation  produced  by  large  aspiration  pressures.  The  non-Newto- 
nian nature  of  the  cell  is  also  suggested  by  its  contents,  for  example  actin  filaments,  microtu- 
bules and  other  suspended  particles,  even  though,  their  role  may  increase  the  cell  viscosity 
only.  It  is  clear  from  experimental  observations,  as  summarized  in  the  previous  section,  that 
the  apparent  viscosity  changes  with  rate  and  extent  of  deformation.  Thus,  the  behavior  of 
neutrophils  cannot  be  described  in  the  same  fashion  as  simple  Newtonian  liquid  drops. 
These  types  of  behaviors  have  been  thought  to  indicate  a  "fading  elastic  memory"  reminis- 
cent of  Maxwell  liquids.  Dong  et  al.  ( 1 988)  use  a  Maxwell  liquid  model  with  a  constant  corti- 


Figure  1 .3.  Illustration  of  a  two-layer  Maxwell  liquid  model  with  constant 
cortical  tension  for  the  leukocyte. 

cal  tension,  as  illustrated  in  Figure  1 .3,  to  analyze  the  behavior  of  neutrophils  flowing  contin- 
uously into  a  pipet.  The  deviatoric  stress  x  is  obtained  in  such  a  model  from: 

^u  +  f^ij  =  2Kj  (1.1) 

where  fX  and  k  are  model  constants  controlling  the  viscous  and  elastic  properties  of  the  Max- 
well liquid.  The  strain  rate  ri  is  obtained  from 

_  I  avj  5Vj 

^iJ  -  2^  ^  (1.2) 

where  v  is  the  fluid  velocity  and  x  represents  spatial  coordinate.  However,  experimental  data 
could  be  correlated  with  their  model  only  if  the  material  properties  of  the  Maxwell  liquid 
were  continuously  varied  as  the  deformation  proceeded.  Hence,  their  model  cannot  be  con- 
sidered to  yield  a  fundamental  description  of  the  leukocyte  in  terms  of  a  Maxwell  fluid. 

Another  attempt  to  obtain  a  non-Newtonian  description  of  the  leukocyte  is  made  by 
Tsai  et  al.  (1993)  through  a  power-law  representation  of  the  cell  apparent  viscosity.  The  ex- 
pression used  by  these  authors  is  of  the  form: 

where  \io  is  the  viscosity  at  a  reference  shear  rate  r\o-  The  constant  b  is  obtained  empirically. 
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The  power-law  model  essentially  tries  to  link  the  shear  rate  dependent  with  the  cyto- 
plasmic viscosity  to  explain  the  shear-thinning  behavior  observed  in  the  micropipet  experi- 
ments. Although  the  power-law  fluid  model  could  be  matched,  for  suitable  initial  condi- 
tions, to  a  fairly  wide  region  of  the  deformation  path  of  leukocytes,  the  rapid  recoil  region 
could  not  be  predicted. 

1.6.  Summary 

Numerical  study  of  leukocyte  modeling  has  relied  on  methods  based  on  semi-analytical 
approaches  (Schmid-Schonbein  et  al.,  198 1 ;  Tran-Son-Tay  et  al.,  199 1)  or  finite-element- 
based  structural  analysis  (Dong  et  al.,  1988;  Dong  and  Skalak,  1992).  It  is  our  opinion  that 
the  full  scope  of  the  numerical  techniques  available  to  the  computational  fluid  dynamics 
community  has  not  been  brought  to  bear  on  the  study  of  leukocyte  modeling. 

It  is  clear  from  the  review  of  current  status  of  leukocyte  modeling  in  the  previous  section 
that  the  existing  viscoelastic  solid,  Newtonian  and  non-Newtonian  liquid  drop  models  fail 
to  explain  the  complex  rheological  behavior  of  leukocytes.  This  is  because  such  models  do 
not  take  into  account  the  internal  structure  of  the  cells.  As  Hochmuth  et  al.  (1993)  point  out 
in  the  conclusion  to  their  review  of  the  state  of  affairs  in  leukocyte  modeling,  it  seems  that 
the  correct  model  of  the  leukocyte  can  only  be  obtained  when  it  is  true  to  the  morphology 
of  the  cell  and  its  contents.  However,  since  the  internal  structure  of  the  cell  is  extremely  com- 
plex, how  can  a  morphologically  consistent  model  be  arrived  at? 

One  strategy,  as  suggested  in  the  three-layer  representation  (Skalak  et  al.,  1990;  Dong 
etal.,  1991;  Hochmuth  etal.,  1993)  shown  in  Figure  1.4,  is  to  simplify  the  situation  by  factor- 
ing in  the  most  important  physical  components,  namely  the  membrane  cortex,  the  cytoplasm 
and  the  nucleus  and  to  temporarily  ignore  the  contributions  of  the  smaller  intracellular  par- 
ticles including  the  microtubules  and  filaments.  This  is  a  good  first  order  approximation 
since  the  cytoplasm  is  mainly  composed  of  water.  The  cytoplasmic  elements  are  assumed 
to  only  produce  a  higher  cytoplasmic  viscosity. 
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Middle  layer  (cytoplasm) 
Figure  1 .4.  schematic  of  a  three-layer  or  compound  drop  model  for  the  leukocyte. 

In  the  three-layer  model  or  compound  drop  under  study  here,  the  outer  layer  is  the  corti- 
cal region  surrounding  the  cytoplasm.  The  cortex  is  comprised  of  an  actin  network.  Its  thick- 
ness is  estimated  to  be  of  the  order  of  0. 1  microns  (Zhelev  and  Hochmuth,  1994)  compared 
to  the  cell  diameter  of  about  8  microns.  Thus,  the  cortical  layer  can  be  considered  to  be  a 
thin  shell  with  constant  tension,  estimated  at  0.024  to  0.038  mN/m  (Evans  and  Yeung,  1989; 
Needham  and  Hochmuth,  1990;  Zhelev  and  Hochmuth,  1994).  The  second  layer,  namely 
cytoplasm,  is  of  a  more  complex  nature  and  under  the  three-layer  viewpoint,  some  scenarios 
have  been  advanced.  It  could  be  a  two-fluid  mixture,  with  a  highly  viscous  liquid  intermixed 
with  a  less  viscous  liquid.  Osmotic  stresses  may  exist  between  these  two  liquids.  When  the 
cell  is  deformed,  the  less  viscous  liquid  would  flow  preferentially.  The  resulting  osmotic 
stresses  would  tend  to  restore  the  equilibrium  state  between  the  two  fluids.  With  regard  to 
the  third  layer,  the  nucleus,  researchers  (Skalak  et  al,  1990;  Dong  et  al.,  1991;  Hochmuth  et 
al.,  1993)  suggest  that  it  can  be  modeled  as  a  highly  viscous,  solid-like  component  which 
retains  its  shape  when  the  cell  is  deformed.  Therefore,  the  leukocyte  is  modeled  as  a  thin 
cortical  shell  with  a  persistent  isotropic  surface  tension  (outer  surface)  surrounding  a  thick 
layer  of  cytoplasm  (shell)  with  a  nucleus  inside  (core  and  inner  interface).  This  representa- 
tion is  also  known  as  a  compound  drop  model. 


11 


The  layout  of  this  thesis  is  as  follows.  In  Chapter  2,  the  details  of  the  salient  features 
of  the  computational  methodology  employed  are  given.  In  Chapter  3,  several  well-defined 
and  challenging  flow  problems  are  tested  to  assess  the  capabilities  of  the  current  numerical 
technique  in  terms  of  accuracy,  robustness  and  flexibilities.  The  large  deformation  and  re- 
covery behavior  of  viscous  compound  drops  are  given  in  Chapter  4.  This  phase  of  the  study 
elucidates  the  hydrodynamics  of  compound  drops  in  relation  to  simple  Newtonian  drops  and 
provides  a  basis  for  the  development  of  a  three-layer  leukocyte  model.  In  Chapter  5,  the 
hydrodynamics  of  the  three-layer  model  and  its  implications  for  the  rheological  behavior 
of  the  leukocytes  are  discussed.  Finally,  Chapter  6  summarizes  the  findings  of  this  study  and 
provides  a  list  of  suggestions  for  future  research. 


CHAPTER  2 
THEORETICAL  ANALYSIS 

2.1.  Background 

The  challenges  posed  by  fluid  flow  problems  involving  moving  boundaries  have  re- 
sulted in  a  variety  of  computational  techniques  directed  toward  their  solution  (Shyy  et  al., 
1996).  The  moving  boundaries  may  represent  material  discontinuities  in  the  flow  field,  and 
their  dynamics  is  coupled  with  that  of  the  flow.  Numerical  methods  applied  to  such  problems 
are  required  to  follow  the  evolution  of  the  often  complex  shapes  assumed  by  the  moving 
boundaries.  To  tackle  these  evolving  interfaces,  as  a  natural  extension  of  adaptive  grid  meth- 
ods for  stationary  complex  geometries,  algorithms  based  on  body-fitted  coordinates  have 
been  applied  to  moving  boundary  problems  (Kang  and  Leal,  1987;  Glimm  et  al.,  1988).  A 
boundary-conforming  grid  arrangement  is  convenient  and  accurate;  the  discontinuity  across 
the  interface  is  always  maintained  and  interfacial  conditions  are  applied  at  the  exact  locations 
without  any  smearing  or  redistribution.  However,  when  the  interface  deformation  becomes 
large,  it  is  difficult  to  adequately  resolve  the  geometrical  complexities  while  maintaining  de- 
sired mesh  control.  Mesh  skewness  and  stretching  may  impact  negatively  on  accuracy  and 
efficiency  of  the  flow  solver.  The  boundary-fitted  grid  method  presents  further  difficulties 
when  interfaces  merge  or  fragment. 

In  simulating  multi-layer  fluid  phenomena  involving  highly  deformed  interface 
shapes,  purely  Eulerian  methods  (Sethian,  1990;  Gunstensen  et  al.,  1991;  Kothe  and  Mjols- 
ness,  1992;  Grunau  et  al.,  1993;  Sussman  et  al.,  1994)  have  come  to  be  widely  used.  In  such 
methods,  a  fixed  Cartesian  grid  is  usually  used,  so  problems  associated  with  grid  generation 
are  circumvented.  Eulerian  methods  based  on  the  Volume  of  Fluid  (Hirt  and  Nichols,  198 1 ), 
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Level-Set  (Sethian,  1990)  or  Lattice  Boltzmann  Equation  (Grunau  et  al.,  1993)  have  been 
developed  to  handle  free-surface  flows,  such  as  wave-breaking,  sloshing  and  splashing  type 
problems,  and  multi-phase  fluid  flows  in  complex  geometries.  In  purely  Eulerian  methods, 
the  details  of  the  interface  shape  are  not  explicitly  needed,  since  the  interface  is  deduced 
based  on  a  computed  field  variable.  The  interface  information  can  be  obtained  from  an  ap- 
propriate computed  field  variable  in  order  to  describe  the  details  of  the  interface  shape  and 
curvature.  An  accurate  estimation  of  the  interface  curvature  is  very  important  when  capillar- 
ity plays  a  significant  role  in  the  interfacial  physics.  In  tracking  highly  distorted  interfaces 
with  good  accuracy,  a  combination  of  the  strengths  of  Eulerian  and  Lagrangian  methods  has 
been  shown  to  be  useful.  This  mixed  Eulerian-Lagrangian  approach  has  been  used  in  recent 
years  for  problems  involving  fluid-fluid  (Unverdi  and  Tryggvason,  1992;  Nobari  et  al., 
1996)  as  well  as  solid-liquid  interactions  (Udaykumar  and  Shyy,  1996;  Juric  and  Tryggva- 
son, 1996).  These  methods  employ  fixed  structured  grids,  and  afford  sifnplicity  and  avail- 
ability of  well-tested  flow  solvers.  The  only  moving  component  is  the  interface  which  is 
a  lower-dimensional  surface  overlying  the  fixed  grid  and  is  explicitly  tracked. 

One  limitation  of  the  methods  based  on  fixed  grids  is  that  since  they  have  traditionally 
employed  Cartesian  grids,  complex  flow  geometries  are  difficult  to  incorporate.  A  complex 
geometry  facility  is  required  for  studying  the  interaction  of  the  moving  boundary  with  flows 
in  arbitrary  domains.  The  representation  of  a  solid  boundary  on  a  fixed  Cartesian  grid  layout 
calls  for  special  treatment,  since  the  control  volumes  through  which  the  irregular  boundary 
passes  become  fragmented  into  solid  and  liquid  regions.  This  situation  is  illustrated  in  Figure 
2.1.  Quirk  (1992)  details  some  of  the  procedures  involved  in  dealing  with  cell  fragments 
or  cut  cells  in  the  framework  of  compressible,  inviscid  flows  around  stationary  obstacles. 
However,  the  calculation  of  fluxes  is  not  presented  in  detail  in  that  work.  Also  in  connection 
with  compressible,  inviscid  flows,  Pember  et  al.  ( 1 995)  use  local  refinement  and  a  redistribu- 
tion procedure  to  enforce  conservation  in  the  vicinity  of  the  solid-fluid  boundary.  Arbi- 
trary-shaped boundaries  running  through  a  Cartesian  mesh  have  also  been  dealt  with  using 
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Figure  2. 1 .  Illustrations  of  (a)  the  situation  when  a  solid  object  in  represented 
on  a  Cartesian  grid,  (b)  a  typical  control  volume  encountered  in  the  vicinity 
of  the  solid-liquid  interface  (Interfacial  cell). 


finite  element  methods  (Young  et  al.,  1991).  In  Zeeuw  and  Powell  (1990)  and  Bayyuk  et 
al.  (1993),  a  Cartesian  grid  is  employed  to  track  the  motion  of  solid  objects  through  an  invis- 
cid  compressible  fluid.  Goldstein  etal.  (1993, 1995)  employ  a  modified  immersed  boundary 
technique  (Peskin,  1977)  to  impose  no-slip  boundary  conditions  over  arbitrary-shaped  solid 
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boundaries.  This  facility  is  useful  for  complex  geometries,  an  important  consideration  for 
the  spectral  methods  in  use  therein.  Turbulent  flows  around  solid  obstacles  have  been  calcu- 
lated using  this  approach  and  the  results  have  been  shown  to  compare  favorably  with  experi- 
ment. The  method,  as  used  in  these  works,  treats  the  solid  boundary  as  a  series  of  steps  con- 
forming with  the  Cartesian  grid  layout,  i.e.,  in  piecewise  constant  fashion.  In  contrast,  a 
piecewise  linear  representation  of  the  free  surface  is  used  by  Miyata  (1986)  for  the  simula- 
tion of  wave  breaking;  the  calculations  are  performed  on  a  Cartesian  grid  using  a  finite-dif- 
ference technique.  In  a  finite  volume  setting,  Udaykumar  and  Shyy  (1995a)  use  a  cut-cell 
technique  to  compute  incompressible,  viscous  flows  in  arbitrary  geometries.  In  that  work, 
it  is  demonstrated  that  by  means  of  a  suitably  defined  consistent  mosaic  of  control  volumes, 
it  is  possible  to  obtain  accurate  solutions  for  flows  on  Cartesian  meshes.  However,  the  stag- 
gered variable  layout  in  Udaykumar  and  Shyy  (1995a)  renders  the  measures  adopted  to  han- 
dle the  internal  boundaries  rather  tedious. 

In  this  work,  we  present  a  fixed  grid  methodology  for  the  simulation  of  flow  through 
complex  geometries  with  fluid-fluid  moving  boundaries,  considerably  simplified  by  the 
choice  of  a  collocated  variable  layout  (Lien  and  Leschziner,  1994).  The  performance  of  the 
flow  solver  in  terms  of  stability  and  accuracy  is  thoroughly  tested  for  various  flow  condi- 
tions. 

2.2.  ELAFINT:  A  Mixed  Eulerian-Lagrangian  Solution  Algorithm 

2.2. 1 .  Characterishcs  of  the  Numerical  Methodology 

There  are  two  distinct  aspects  to  the  computational  framework,  called  ELAFINT,  which 
stands  for  Eulerian-Lagrangian  Algorithm  for  INterface  Tracking,  reported  in  Udaykumar 
et  al.  (1995a,  1996).  These  are: 

(1)  A  pressure-based  flow  solver  to  simulate  the  incompressible  fluid  flow  over  a  wide 
range  of  Reynolds  numbers. 
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(2)  A  procedure  to  deal  with  the  interaction  of  fixed  arbitrary  geometries  and  evolving 
moving  boundaries  with  the  flow. 

The  first  aspect  is  straightforward  and  approached  in  a  conventional  manner  using  the 
robust  pressure-based  solution  procedure  detailed  by  Patankar  (1980)  and  in  the  context  of 
more  modern  developments  by  Shyy  (1994).  The  scheme  to  track  the  moving  fronts  and 
describe  complex  geometries  is  summarized  below,  following  which  the  flow  solver  and  its 
discretization  procedure  are  described. 

2.2.2.  Interface  Tracking  Procedure 

The  interfaces  are  described  by  means  of  markers  indexed  sequentially  and  separated 
into  objects.  There  is  no  restriction  on  the  number  of  objects  or  on  the  end  conditions  that 
can  be  imposed  on  these  curves.  As  shown  in  Figure  2.2,  the  objects  can  be  open  or  closed. 
Furthermore,  they  can  be  fixed  or  moving.  They  can  enclose  different  phases  or  materials 
with  different  properties.  Once  the  interface  description  is  obtained  in  terms  of  markers  and 
their  positions,  the  shape  characteristics  such  as  normals  and  curvatures  can  be  extracted. 

The  convention  adopted  by  the  algorithm  is  to  consider  the  normal  to  the  interface  as 
pointing  from  phase  II  to  phase  I,  as  shown  in  Figure  2.2,  such  that  when  one  traverses  the 
interface  along  the  marker  sequence,  phase  II  lies  to  the  right.  In  order  to  handle  multiple- 


j+1 
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n 


PHASE  I 


PHASE  II 


Figure  2.2.  Illustration  of  objects,  markers  and  normal  convention 
in  the  description  of  interfaces. 
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valued  interfaces,  the  interface  tracking  algorithm  employs  piecewise  polynomial  fits  to 
compute  the  necessary  shape  derivatives  at  each  marker  location.  Facility  exists  to  construct 
quadratic  or  cubic-spline  functions  Xi(s)  andyi(s).  The  polynomial  fits  are  obtained  by  para- 
metrizing the  curve  as  a  function  of  arc  length  s  so  that,  for  the  ith  marker 

Xi  =  Xi(s),     y,.  =  yiis)  (2.1) 

Then,  the  normal  and  curvature  are  obtained  by  taking  appropriate  derivatives  of  these  piece- 
wise  polynomials.  Thus,  the  components  of  the  normal  are, 

^  '-^   ^  


+  ^^2)1/2  '        y  (x,2-Ky/)l/2  (2.2) 

The  curvature  x  is,  for  the  two-dimensional  case: 

•  n  =   — TTT 

O'.^  +  Xs^)y^  (2.3) 

and  the  axisymmetric  case: 

   Xs  |_     ysXsS  ~ 

"  y(ys^  +  Xs^)y^  +  x,2)3/2  (2.4) 

Interfacial  markers  are  advected  to  their  new  positions  by  a  Lagrangian  translation.  Thus, 

=       +  At  Uf  (2.5) 

y/+^  =  y^l^+AtVf  (2.6) 

where  the  subscript  £  indicates  the  interfacial  marker  number  and  the  superscript  k  implies 
the  time  level.  It  is  demonstrated  (Udaykumar  and  Shyy,  1995a,  b;  Shyy  et  al.,  1996)  that 
the  interface  tracking  procedure  can  handle  extreme  deformation,  including  topological 
changes  such  as  repeated  mergers  and  breakups.  It  is  also  established  that  the  interface  track- 
ing facility  is  versatile  and  accurate. 

Once  the  objects  are  described  on  the  computational  domain,  the  interaction  with  the 
flow  field  needs  to  be  established.  This  is  mainly  done  by  providing  the  following  informa- 
tion: 
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(1)  The  cell  in  which  each  interface  marker  lies. 

(2)  The  markers  contained  within  each  interfacial  cell  (i.e.,  control  volume  through  which 
the  interface  passes). 

Since  a  Cartesian  grid  is  employed,  establishing  these  relationships  is  a  straightforward 
task.  Based  on  the  information  mentioned  above,  the  flow  solver  can  incorporate  either 
changes  in  the  control  volume  definitions  (as  in  the  cut-cell  approach)  or  appropriate  source 
term  contributions  (as  in  the  immersed  boundary  technique).  These  two  approaches  enable 
the  effects  of  the  internal  stationary  as  well  as  moving  boundaries  to  be  communicated  with 
the  flow  field.  In  turn,  based  on  the  flow  field  developed,  the  motion  of  the  interface  can 
be  determined. 

In  Udaykumar  and  Shyy  (1995b),  moving  as  well  as  stationary  solid-liquid  boundaries 
are  handled  by  the  cut-cell  technique.  In  this  work,  the  cut-cell  technique  is  retained  for  the 
stationary  solid-liquid  boundaries,  while  the  immersed  boundary  technique  (Peskin,  1977; 
Fauci  and  Peskin,  1988;  Unverdi  and  Tryggvason,  1992)  is  used  to  treat  moving  fluid-fluid 
boundaries.  This,  along  with  the  use  of  a  collocated  variable  arrangement,  enhanced  the  ro- 
bustness and  simplicity  of  the  algorithm.  Extensive  testing  of  the  algorithm  has  been  per- 
formed for  a  wide  variety  of  flows  and  the  results  thereof  are  presented  elsewhere  (Shyy  et 
al.,  1996;  Udaykumar  et  al.,  1996). 

2.2.3.  Formulation  of  the  Pressure-Based  Flow  Solver 

Consider  the  two-dimensional  planar  or  axisymmetric,  incompressible  Navier-Stokes 
flow.  The  governing  equations  in  conservation  forms  are: 
Continuity  equation: 

^  +  ^(PM)  +  -L^(y»  =  0 
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Momentum  equations: 


d{gu) 
dt 


dx  ^  dx^  dx' ^  y'"  dy^y  ^df  ^  ^' 


u 


(2.8) 


d{Qv)  _^  _d_ 
dt  dx 


dy     ^  ^2  +  dx^dx^  ^  y^dy^  ^  dy^  ^ 


(2.9) 


Energy  equation: 


dt  dx 


(2.10) 


where  m=0  for  planar  and  1  for  axisymmetric  case.  In  the  above,  g  is  the  fluid  density,  u 
and  V  are  the  fluid  velocity  components,  T  is  the  temperature,  p  is  the  pressure  and  fi  is  the 
coefficient  of  viscosity  and  5„  and  5^  are  the  source  terms  in  x-  and  y-momentum  equations 
which  includes,  for  example,  the  buoyancy  term  under  the  Boussinesq  approximation  and 
the  surface  tension  contributions  discussed  later. 

The  boundary  conditions  to  be  applied  on  the  interface  for  cases  involving  no  mass  ex- 
change between  phases  are  the  continuity  condition  and  normal  stress  balance: 


where  the  subscripts  1  and  2  represent  the  two  adjoining  fluids  and  /  represents  the  interface, 
V„  is  the  normal  velocity  component,  y  is  the  surface  tension  and  x  is  the  curvature  of  the 
interface,  with  the  normal  pointing  from  fluid  2  into  fluid  1,  as  shown  in  Figure  2.3. 


(V„),  =  iVn)2  =  (Vn), 


(2.11) 


(2.12) 
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normal  to  interface 


Underlying  fixed 
Cartesian  grid 


Liquid  2 


normal  to  solid  surface 


Figure  2.3.  Ilustration  of  the  interaction  of  the  solid-liquid  and 
liquid-liquid  boundaries  with  Cartesian  grid.  The  sohd-liquid  is  treated 
using  cut-cell  technique,  while  the  liquid-liquid  interface  is  modeled 
with  the  immersed  boundary  technique. 

2.2.4.  Discretization  Procedure  of  the  Flow  Field  Equations 

The  discretization  of  the  governing  equations  is  carried  out  using  a  control  volume  for- 
mulation, which  with  particular  reference  to  two-dimensional  planar  flow.  Similar  treat- 
ment applies  for  the  axisymmetric  case.  Consider  the  conservation  law  for  the  variable  0, 
defined  to  be  (a)  1  for  the  continuity  equation  (b)  u  and  v  for  the  x-  and  y-  momentum  equa- 
tions and  (c)  specific  enthalpy  (h=CpT)  for  the  energy  equation.  In  the  general  case,  we  have 


The  control  volumes  in  the  vicinity  of  a  solid  boundary  passing  through  the  grid  are 
irregularly  shaped,  and  in  the  general  case  can  assume  a  five-sided  shape  as  shown  in  Figure 
2.1(b).  In  order  to  evaluate  the  fluxes  through  the  cell  faces,  we  integrate  Eq.  (2.13)  over 


d(Q(p) 
dt 


+  V  •  (^M0)  =  (V  •  rV(^)  +  s 


(2.13) 
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the  control  volume  and  employ  the  divergence  theorem  for  two-dimensional  problems  to 
obtain 


dt 


dA  + 


-0)  •  ndl   =   ()(rV0)  •  ndl  + 


I 


SdA 


I  I  A  (2.14) 

where  the  second  term  on  the  left  hand  side  is  now  the  line  integral  of  the  outward  normal 
convective  fluxes  and  the  first  term  on  the  right  side  is  the  line  integral  of  the  outward  normal 
diffusion  fluxes  through  the  faces  of  the  control  volume. 

We  now  proceed  to  discretize  each  of  the  terms  in  Eq.  (2.14)  for  the  control  volume 
shown.  Thus 


1 


k+\^k+\    _  Qkfkk 


dt 


(2.15) 


where  the  superscript  k,  k+1  indicate  the  time  levels  and  dt  is  the  time  step  size.  A^v  is  the 
area  of  the  irregular  control  volume.  The  time  derivative  is  discretized  in  the  backward  Euler 
form  leading  to  an  implicit-in-time  solution  procedure  which  relaxes  the  time  step 
constraints  in  comparison  with  explicit  schemes.  Next, 


5 


(){qu<P)  •  ndl    =    y  (QfU^f) 


k+\ 


dip 


€=1 


/  -   •  (2.16) 

is  the  summation  of  the  convective  normal  effluxes  through  each  cell  face  of  the  control  vol- 
ume. The  superscript  k-^l  indicates  the  implicit  nature  of  the  scheme.  Following  Shyy  et 
al.  (1996)  and  Udaykumar  et  al.  (1996),  the  diffusion  fluxes  are  computed  as 

5 


Lni4>)-ndi  =  Y.'<^((^h)'^'di 
i 


(2.17) 
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Let  us  denote  the  source  term  as  follows, 


1 


SdA    =  S 

A  (2.18) 
Here  S  contains  contributions  from  the  pressure  and  body  force  terms.  Substituting  the 
above  Eqs.  (2.15-2.17)  in  Eq.  (2.14)  one  obtains  the  discretized  form  as 

t=\  e=\  (2.19) 

In  particular  the  discrete  form  of  the  continuity  equation  can  be  written  as 

(e)*+'-p*.  )  5 

"    .       "    Ac,    +     Y.{QUn)^dl^    =  0 

(  =  1  (2.20) 
Now  multiplying  Eq.  (2.20)  by  the  value  <pij  and  subtracting  from  Eq.  (2.19)  gives 

-^^^^^  ^A,,  +  2^  Q,u,i<p,  -  <t>i/^'dl,  =  ^(r,(£),)*+ij/,  +  S 

e  =  \  €=1  (2.21) 

In  the  notation  of  the  pressure-based  methodology  (Patankar,  1980;  Shyy,  1994),  we  cast 
the  final  discretized  form  as 

where 

dd) 

b    =   S  +  r l-^)dl,  -  Qiu„,dl^<p,  -  (pp) 

Thus,  a  five-point  stencil  is  adopted,  and  the  value  of  the  variable  (p  at  point  P  is  dependent 
on  its  four  adjoining  neighbors.  The  form  of  the  coefficients  ciewns  depends  on  the  convec- 
tion scheme  employed.  In  this  work,  we  employ  the  second-order  central  difference  scheme 
for  the  convection  term.  Diffusion  and  pressure  terms  are  also  discretized  using  central  dif- 
ferences. When  cast  in  the  matrix  form  the  above  equation  reads 

[C]  m    =   [B]  (2.24) 
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where  [C]  is  a  pentadiagonal  coefficient  matrix,  0  is  the  solution  vector  and  B  is  a  source 
vector.  The  procedure  usually  adopted  is  to  solve  the  system  of  equations  in  iterative  fashion 
employing  a  line  SOR  method  which  calls  for  a  tridiagonal  matrix  solution  procedure.  A 
pressure-correction  procedure  provides  the  pressure  field  consistent  with  mass  conserva- 
tion. In  the  pressure-based  methodology,  it  is  common  practice  to  use  a  staggered  layout 
of  variables.  However,  here  we  use  a  collocated  arrangement  of  the  dependent  variables. 
This  calls  for  special  measures,  called  momentum  interpolation  (Rhie  and  Chow,  1983;  Lien 
and  Leschziner,  1 994),  when  estimating  cell  face  velocities.  The  details  of  momentum  inter- 
polation will  be  discussed  next. 

2.3.  Momentum  Interpolation 

One  of  the  key  issues  in  using  a  pressure-based  solver  is  the  grid  arrangement.  The 
staggered  grid  layout  has  been  preferred  for  incompressible  flow  computations  for  many  de- 
cades due  to  its  many  advantages  (Patankar,  1980;  Shyy,  1994).  However,  the  collocated 
arrangement  has  desirable  attributes  in  terms  of  simplicity  and  economy,  has  come  into  com- 
mon use.  For  example,  in  the  present  work,  in  dealing  with  the  internal  boundaries,  whether 
by  the  cut-cell  or  the  immersed  boundary  technique,  the  use  of  a  collocated  variable  layout 
leads  to  the  treatment  of  a  single  set  of  control  volumes  for  all  the  dependent  variables  (name- 
ly, u,  V,  p  and  T).  In  the  case  of  a  staggered  system,  three  sets  of  control  volumes,  for  u,  v 
and  the  scalars  need  to  be  considered  and  the  data  structures  associated  with  the  cut-cell  as 
well  as  immersed  boundary  technique  need  to  be  computed  and  stored.  A  particular  instance 
is  the  identification  of  the  phase  of  the  control  volume. 

When  the  staggered  layout  is  used  the  phase  in  which  the  control  points  for  the  u,  v  and 
scalar  variables  lie  needs  to  be  determined  separately,  while  in  the  collocated  layout  only  one 
set  of  control  volumes  needs  to  be  handled.  This  leads  to  a  considerable  savings  both  in  terms 
of  the  programming  effort  and  also  in  the  computational  effort  expended  in  obtaining  the 


24 


information  regarding  the  interaction  of  interfaces  with  the  flow  field.  These  two  variable 
arrangements  are  illustrated  in  Figure  2.4. 

One  potential  difficulty  of  the  collocated  grid  is  the  occurrence  of  unrealistic  oscilla- 
tions in  the  pressure  and  velocity  fields,  leading  to  numerical  instabilities  for  certain  flow 
problems  (Patankar,  1980).  This  is  due  to  the  fact  that  the  collocated  grid  arrangement  leads 
to  decoupling  of  pressure  and  velocity  fields  when  the  central  difference  operator  is 
employed  in  discretizing  gradients  in  the  flow  field.  However,  with  the  incorporation  of  mo- 
mentum interpolation,  introduced  by  Rhie  and  Chow  (1983),  it  is  possible  to  prevent  the  on- 
set of  checkboard  phenomena  resulting  from  the  decoupling  between  the  pressure  and  veloc- 
ity fields.  Lien  and  Leschziner  (1994)  have  applied  this  measure  to  obtain  solutions  for  a 
wide  range  of  flow  problems,  including  shocked  flows.  In  our  implementation,  we  use  the 
version  of  "momentum  interpolation",  recommended  by  Lien  and  Leschziner  ( 1 994)  to  eval- 
uate the  velocities  at  the  faces  of  the  control  volume. 

Essentially,  momentum  interpolation  introduces  a  third-order  dissipation  in  estimating 
the  cell  face  velocity.  This  can  be  seen  as  follows. 

The  discretized  momentum  equations  for  up  and  ue  based  on  finite  volume  formulation, 
presented  for  a  general  transport  variable  (p  in  Eq.  (2.13)  are  given  in  the  one-dimensional 
case  by 

ApUp  =     ^    a„^M„^  +  Sp-  ipe-  Pw)p 

nb  =  E,W  (2.25) 

^£«£  =      X     ^"bl^nb  +  ^E-  (Pe~  Pw)e 

nb  =  EE,P  (2.26) 

where  subscript  nb  implies  neighboring  points  in  the  stencil,  Sp  and  Se  represent  the  source 
terms,  Ap  =    ^         and  Ap-  -     ^    a^^.  The  relative  positions  of  the  points  sub- 

nb  =  E,W  nb  =  EE,P 

scripted  P,  E,  W,  EE  and  WW  are  shown  in  the  Figure  2.5. 
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Figure  2.4.  Schematics  of  (a)  staggered  grid  layout,  (b)  collocated  grid  layout. 
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Figure  2.5.  Schematic  of  1-D  collocated  grid  layout. 


Equation  (2.25)  and  Eq.  (2.26)  can  be  rewritten  as 

Hp 


-  4-(Pe  -  Pw)p 


(2.27) 


1  , 

"e  =  7  -J-^e  ~  P^>E 


(2.28) 


where =    ^   ^nb'^nb     SpSindHg  =     ^    ^nb^nb     ^E-  In  performing  the  inter- 
na =£,w  nb  =  EE,P 

polation,  it  is  assumed  that 

He^\lHp^HE\ 

Ag     2\Ap     Ai;j  ^2.29) 

Then,  substituting  Eq.  (2.27)  and  Eq.  (2.28)  into  Eq.  (2.39),  the  interpolation  for  the  face 
velocity  Ug  can  be  expressed  as 


or  alternatively. 


1 


1 


~  ~  Pe)p  +  «£  ~  T~^^  ~  Pe)E 


+  ~  Pe) 


(2.30) 


Ue  =  ^(Mp  +  Mg)  + 


-^{pp  -  Pe)-\ 


'^(Pw  ~  Pe)p  +  '^^^  ~  Peh 


(2.31) 


In  Eq.  (2.31)  the  face  velocity  of  control  volume  is  directly  linked  to  the  two  adjacent 
cell-center  pressures  through  the  term  I  above.  This  ensures  that  strong  coupling  between 
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the  pressure  and  velocity  fields  is  retained.  The  second  term,  indicated  as  II,  is  the  pressure 
smoothing  term,  which  eliminates  the  non-physical  wavy  pressure  or  velocity  fields.  In  fact, 
it  can  be  shown  that  Eq.  (2.3 1 )  can  be  recast,  following  some  manipulation  (Lien  and  Lesch- 
ziner  1994),  as 


1  /av' 

Ue  =  -j^iUp  +  Mg)  +  I 


lO(zlx3) 


(2.32) 


This  brings  to  light  the  fact  that  momentum  interpolation  adds  a  third-order  term  to  the 
cell  face  velocity,  relative  to  straightforward  averaging  of  adjacent  cell-center  velocities. 
Thus  spurious  pressure  oscillations  are  eliminated.  Typically,  iterative  procedures  are  used 
to  solve  the  discretized  momentum  equations,  and  under-relaxation  can  be  implemented  in 
Eq.  (2.27)  in  the  form 


Up  =  a  Up  +  {\  —  a)  Up 


(2.33) 


where  a  is  the  under-relaxation  factor,  up  is  the  current  iterative  value,  Mpis  the  value  of  up 
at  the  previous  iteration,  while  u*p  is  the  value  obtained  after  relaxation.  Similarly, 


Ue  =  a  Up  +  {\  —  a)  ue 


Therefore,  the  cell  face  velocity,  after  including  the  relaxation  factor  is 


(2.34) 


Ue  =  a  Ug  +  {\  —  a)  Ue 


=  a- 


^(up  +  Ue)  +  j-^P  ~  Pe) 


■J-iPw  -  Pe)p  +  4-iPw-  Pe)E 
/\p  Zip 


+  (I  —  a)  Ug 


(2.35) 


The  momentum  interpolation  as  obtained  from  Eq.  (2.31)  yields  flow  fields  indepen- 
dent of  the  relaxation  factor,  thus  eliminating  the  relaxation  factor  dependency  detected  in 
early  versions  of  momentum  interpolation  (Majumdar,  1988).  From  Eqs.  (2.33-2.34),  Up 
and  u E  can  be  rewritten  as 
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♦* 

Up 

Up  = 


1  -  a 
a 


Up 


(2.36) 


** 
Ue_ 

a 


1  -  a 

a 


Ue 


(2.37) 


Thus,  Eq.  (2.35)  becomes 


uT  =  7;(u*p*  +  u*e)  +  (1  -  a) 


Ue  -  ^(Mp  +  Ue) 


+  a 


■J-iPw  -  Pe)p  +  -7-iPw  -  Pe)E 


(2.38) 


Equation  (2.38)  is  the  form  of  the  momentum  interpolation  suitable  for  steady  flow 
problems.  The  extension  of  Eq.  (2.38)  for  unsteady  flow  problems  will  include  the  time-de- 
pendent term  {QAxu°/At)  where  superscript  o  indicates  the  previous  time  value.  The  ap- 
propriate form  is 


**      1  M 


Ue    =  2I  J""^   +  Ye^E  I  +  (1  -  a) 


"^"2  U:"^  a7"^ 


Zir  ]  2 


^P  -  Pe)  -  9  [(Pw  -  P^P  +  (Pw  -  Pe)E\ 


Ay 


(2.39) 


Although  Eq.  (2.39)  applies  to  a  1-D,  uniform  grid  situation,  it  can  be  easily  extended 
to  higher  dimensions  and  non-uniform  grids. 


2.4.  Interactions  between  Interfacial  Physics  and  Flow  Field  Variables 


2.4.1.  Cut-Cell  Procedures  for  Solid-Liquid  Boundaries 

In  accounting  for  the  presence  of  the  internal  solid  boundaries  which  define  complex 
geometries  on  Cartesian  grids,  the  cut-cell  approach  developed  by  Udaykumar  and  Shyy 
(1995a)  is  applied.  In  this  approach,  when  the  interface  passes  through  the  grid,  the  scenario 
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shown  in  Figure  2. 1  is  encountered.  The  control  volume  designated  by  P  is  shown  to  be  cut 
by  the  solid  surface  passing  through  the  grid.  The  cut-cell  procedure  rearranges  the  control 
volumes  in  the  vicinity  of  the  interface  to  form  a  fully  consistent  mosaic  of  cells.  This  is 
obtained  by  redefining  the  interfacial  cells  (i.e.,  the  control  volumes  in  the  vicinity  of  the 
interface),  so  that  the  control  volume  fragments  lying  in  each  phase  are  absorbed  into  ap- 
propriate partner  cells  (Udaykumar  and  Shyy,  1996). 


Figure  2.6.  Illustration  of  the  situation  resulting  when  a  complex  solid 
geometry  is  treated  on  a  fixed  Cartesian  grid.  The  cut  interfacial  cells  are 
reassembled  to  provide  a  closed,  consistent  mosaic  to  maintain  flux 
conservation.  The  assembly  of  cells  is  performed  by  absorbing  fragments 
of  the  control  volumes  into  appropriate  cells  to  reconfigure  the  cells  in  the 
vicinity  of  the  boundary  as  shown  by  the  dotted  lines.  Also  shown  in  the 
normal  probe  for  obtaining  gradients  at  the  interface. 

An  example  is  shown  in  Figure  2.6,  where  the  dotted  lines  indicate  the  newly  defined 
control  volumes.  The  fluxes  through  these  redefined  control  volume  faces  are  thereby  eva- 
luated in  such  a  way  as  to  maintain  conservation  and  consistency  of  the  fluxes  across  each 
control  volume  face.  For  example,  the  flux  of  a  conserved  variable  is  computed  for  a  cell 
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that  is  abutted  by  two  neighbors  to  the  left  (as  in  the  case  of  the  hatched  cell  in  Figure  2.6), 
such  that  the  flux  through  the  west  face  of  the  hatched  cell  is  the  sum  of  the  fluxes  through 
the  east  faces  of  the  abutting  cells. 

This  procedure  obviates  any  smearing  of  information  at  the  interface  since  the  boundary 
condition  is  imposed  at  the  exact  location  of  the  internal  boundary  in  each  control  volume. 
This  is  in  contrast  to  the  means  employed  by  Goldstein  et  al.  (1993)  in  imposing  the  no-slip 
condition.  In  that  work,  the  no-slip  boundary  condition  is  imposed  via  an  external  concen- 
trated force  field  which  is  computed  based  on  the  developing  flow  field  in  the  vicinity  of  the 
interface.  Thus,  a  feedback  control  approach  is  used  to  compute  the  shear  stresses  necessary 
to  obtain  a  no-slip  condition.  Information  regarding  the  no-sHp  surface  is  redistributed  over 
the  finite  bound  of  the  d —function  used  in  the  immersed  boundary  representation  of  the  solid 
surface. 

In  this  work,  for  the  evaluation  of  the  source  term  in  Eq.  (2.23),  the  diffusion  flux  from 
the  interface  is  obtained  by  describing  the  field  (p(x,y)  for  any  transport  variable  such  as  ve- 
locity or  temperature.  This  is  done  in  each  phase  separately  around  the  interface  by  means 
of  a  biquadratic  function  and  the  gradient  normal  to  the  interface  required  in  obtaining  the 
diffusion  flux  is  obtained  by  extending  a  normal  probe  from  the  interface  into  each  phase. 
The  length  of  the  normal  is  one  grid  spacing  h.  This  situation  is  illustrated  for  the  liquid  phase 
in  Figure  2.6.  The  value  <p[+dn  at  the  end  of  the  probe  is  obtained  from  the  biquadratic  func- 
tion representation  in  the  required  phase  as  4>(x,y)  =  ax^+by^+cxy+dx+ey+f.  The  coeffi- 
cients are  obtained  by  inverting  the  6x6  matrix  obtained  by  the  known  values  of  (p  at  six  suit- 
ably chosen  points  in  the  required  phase  around  the  end  of  the  probe.  The  inversion  of  the 
matrix  is  performed  by  LU  decomposition.  Thus,  the  required  normal  gradient  of  the  field 
<p(x,y)  is  computed  as  follows: 


(2.40) 
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where  dn  is  the  normal  probe  length.  The  pressure  at  the  solid  boundary  is  obtained  by  com- 
puting the  pressure  in  a  similar  manner  to  the  above  at  the  reference  point  and  extending  the 

same  to  the  interface,  which  implies  that  ^  =  0  at  the  solid  boundaries.  In  the  case  of  the 

stationary  solid  boundaries  in  this  work,  the  convective  interfacial  fluxes  do  not  contribute 
to  the  momentum.  However,  for  moving  solid  boundaries,  the  treatment  is  given  in  Udayku- 
mar  et  al.  (1996).  The  procedures  for  the  treatment  of  the  internal  solid  boundaries  lead  to 
stable,  accurate  computations  of  flows  through  complex  geometries. 

2.4.2.  Immersed  Boundary  Technique  for  Fluid-Fluid  Interfaces 

When  the  fluid-fluid  interface  (i.e.,  the  moving  boundary)  lies  on  a  fixed  grid,  some 
means  of  communication  between  the  interface  and  the  Eulerian  framework  for  the  field 
equations  needs  to  be  devised.  One  such  technique  is  the  cut-cell  method  described  above 
and  in  more  detail  in  Udaykumar  et  al.  (1996).  In  this  work,  we  adopt  the  immersed  bound- 
ary technique  originated  by  Peskin  (1977),  and  employed  extensively  by  Tryggvason  and 
coworkers  (Unverdi  and  Tryggvason,  1992;  Juric  and  Tryggvason,  1996;Nobarietal.,  1996) 
for  a  variety  of  two-phase  flow  problems. 

In  the  immersed  boundary  technique,  the  presence  of  the  interface  and  its  effect  on  the 
flow  field  is  established  by  means  of  source  terms  in  the  governing  equations,  i.e.,  the  in- 
formation from  the  Lagarangian  interface  component  is  redistributed  to  the  Eulerian  frame- 
work. For  the  particular  case  of  two-fluid  interfaces  under  constant  interfacial  tension,  the 
source  term  represents  the  normal  stress  contribution  from  the  interface.  The  momentum 
jump  conditions,  given  by  Eqs.  (2. 1 1-2. 12)  at  the  interface,  is  incorporated  in  the  x-  and  y- 
momentum  equations.  Specifically,  the  contributions  to  the  normal  stress  balance  from  the 
capillarity  term  in  Eq.  (2. 12)  is  redistributed  to  the  surrounding  grid  points;  it  is  supplied  to 
the  source  term  Su  and  in  Eqs.  (2.8-2.9).  This  is  done  with  the  aid  of  a  6  function  of  bound 
4h  as  illustrated  in  Figure  2.7  (where  h  is  the  grid  spacing).  This  procedure  leads  to  the  con- 
servative redistribution  of  the  source  terms  to  the  grid,  due  to  the  properties  of  the  delta  func- 
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tion.  The  discrete  delta  function  at  each  point  in  the  flow  field  assumes  the  form  originally 
proposed  by  Peskin  (1977): 


where  x  and  (x)  ^  are  the  locations  of  points  in  the  flow  field  and  on  the  interface,  respective- 
ly, €  represents  the  interfacial  marker  index.  R  is  set  equal  to  2h  in  our  computations. 

Each  marker  contributes  to  the  momentum  source  term  in  its  proximal  control  points, 
i.e.,  those  lying  in  the  circle  of  radius  2h  centered  at  the  marker.  The  control  points  experi- 
encing the  influence  of  the  capillary  forces  are  marked  in  Figure  2.8  by  closed  dots.  The 
force  contribution  from  each  interfacial  marker  to  the  momentum  equations  is  at  any  point 
P  given  by 
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Figure  2.8.  Illustration  of  the  form  of  the  delta  function  in  the  vicinity  of  the 
interface,  h  is  the  grid  spacing.  The  bound  of  the  delta  function  is  shown. 


fpt  -  y^e^e^^e  (2.42) 

where  zj^f  is  the  arc  length  corresponding  to  the  marker.  Thus,  the  source  term  in  the 
momentum  equations  at  a  given  point  in  the  flow  field  depends  on  the  proximity  of  that  point 
to  the  fluid-fluid  interface. 

The  magnitude  of  the  interfacial  forces  experienced  is  dictated  by  the  delta  function  and 
is  the  sum  of  the  contributions  of  all  the  interfacial  markers  lying  within  the  distance  of  2h 
from  the  point.  Therefore,  the  normal  stress  contribution  to  the  momentum  equation  is  given 
by  the  expression, 

(  (2.43) 


34 


This  procedure  of  redistribution  of  interfacial  stresses  as  sources  to  the  field  equations 
leads  to  the  smoothing  of  the  discontinuity  at  the  interface.  The  stability  and  robustness  of 
this  method  is  enhanced  by  such  smoothing. 

The  continuity  condition  at  the  interface  for  material  interfaces,  Eq.  (2. 1 1),  is  enforced 
by  extracting  the  normal  velocity  of  the  interfacial  markers  from  the  grid  information  by 
means  of  a  delta  function.  The  d  function  form  is 

J  (2.44) 

where  the  index  j  corresponds  to  the  proximal  (in  the  region  given  by  R=2h)  grid  points.  The 
interface  is  advected  with  this  velocity. 

Another  feature  that  needs  to  be  mentioned  is  the  handling  of  property  jumps  across  the 
interface.  In  other  approaches  using  a  fixed  grid  (Kothe  et  al.,  1992;  Sussman  et  al.,  1994) 
some  field  variable  is  employed  to  smooth  the  property  jumps  across  the  interface.  For  ex- 
ample, a  Heaviside  function  can  be  employed  such  that: 


3G(0)a  = 


1,  <p  <  -  a 

0,  <p  >  +  a 

4    a  +  ^^1"—  otherwise 


(2.45) 

where  a  is  the  bound  region  around  the  interface  which  is  equal  to  2h  in  our  computations. 
Usually  this  bound  is  taken  to  be  proportional  to  the  grid  size.  Based  on  the  Heaviside  func- 
tion, the  material  properties  are  assigned  as  : 

^  =/3,  +O32-/3,)aG(0)«  (2.46) 

where  ^  is  any  material  property  such  as  viscosity  and  density.  The  subscripts  2  and  1  indi- 
cate the  two  adjoining  phases,  in  line  with  our  convention  for  the  normal. 

In  the  RIPPLE  approach  (Kothe  et  al.,  1992),  which  is  based  on  a  VOF  representation 
of  the  phase  field,  0  is  taken  as  a  color  function  dependent  on  the  fluid  fraction  field,  while 
in  the  level-set  algorithm  the  ((>  value  corresponds  to  a  distance  function.  In  these  methods 
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the  field  variable  0  is  an  essential  part  of  the  interface  tracking  procedure.  In  the  immersed 
boundary  technique,  Unverdi  and  Tryggvason  (1992)  employ  an  indicator  function  which 
is  obtained  by  solving  a  Poisson  equation  for  the  entire  domain.  The  indicator  function  has 
the  characteristics  of  the  Heaviside  function  but  has  to  be  computed  every  time  the  interface 
moves.  Thus,  the  calculation  of  an  indicator  function  involves  the  solution  of  an  additional 
field  equation  and  does  not  fully  exploit  the  available  information  from  the  explicitly  tracked 
interface.  In  our  algorithm  we  have  exploited  the  Heaviside  representation  to  achieve  the 
property  smoothing,  based  on  explicit  knowledge  of  the  interface  and  the  information  indi- 
cated in  Section  2.2.1. 

Once  the  marker  locations  with  respect  to  the  grid  are  known,  the  property  jumps  can 
be  assigned  based  on  the  Heaviside  function.  In  the  present  case,  the  Heaviside  function 
takes  the  form: 


1,  sgn(x  —  (x)f)\x  —  {x)f\  <  —  a 

0,  sgn(x  —  (x)f)\x  —  (x)^\  >  +  a 


(2.47) 

where  x,  represents  the  coordinate  of  the  control  point  where  the  Heaviside  function  is  re- 
quired. The  smoothing  of  the  material  properties  is  performed  with  Eq.  (2.47)  above. 

The  material  properties  in  the  bulk  of  the  phases  are  assigned  first.  Thereafter  the  transi- 
tion region  between  any  two  phases  is  smoothed  by  defining  the  Heaviside  function  locally 
in  this  region.  To  accomplish  this,  one  traverses  along  the  interface  and  modifies  the  value 
of  the  material  property  according  to  Eq.  (2.47),  in  the  cells  adjacent  to  the  interface  and  ly- 
ing within  the  bound  a  of  the  interface.  These  operations  are  confined  to  the  region  close 
to  the  interface  only,  hence  reducing  the  computational  burden.  This  procedure  for  smooth- 
ing material  discontinuities  holds  when  there  are  multiple  interfaces,  including  cases  in 
which  interfaces  may  enclose  other  interfaces. 
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It  is  noted  that  computational  difficulties  such  as  spurious  oscillations  and  numerical 
instabilities  are  often  encountered  when  the  property  gradient  across  the  interface  is  steep. 
By  applying  the  Heaviside  function  to  smooth  the  discontinuities  across  the  interface,  the 
choice  of  bound  region  a  needs  to  be  made  with  care.  A  large  a  can  cause  the  solution  to 
be  overly  smeared  and  to  substantially  lose  physical  fidelity  while  a  small  a  will  be  insuffi- 
cient to  stabilize  the  numerical  solution.  Another  issue  is  the  selection  of  the  smoothing 
function.  Instead  of  choosing  the  smoothing  function  which  is  linearly  proportional  to  the 
mesh  size,  a  higher  order  smoothing  function  can  be  used  (Grunau  et  al.,  1993). 

2.5.  Mass  Conservation  in  Individual  Phases 

In  the  scheme  presented  above,  the  interaction  of  the  interface  tracking  algorithm  as 
provided  by  the  immersed  boundary  technique,  with  the  flow  solver  does  not  enforce  con- 
servation of  mass  in  the  phase  enclosed  by  the  interface.  This  is  because  the  divergence-free 
velocity  field  is  enforced  in  each  control  volume  by  the  flow  solver,  the  velocity  employed 
to  evolve  the  interface  is  interpolated  from  such  a  velocity  field.  This  interpolated  velocity 
field  may  violate  the  divergence-free  condition,  resulting  in  a  mass  leak  from  the  phase  en- 
closed by  the  interface. 

Peskin  and  Printz  ( 1 993)  suggest  a  method  based  on  the  modification  of  the  divergence 
operator  in  order  to  enforce  mass  conservation.  In  this  work,  we  choose  to  achieve  the  same 
result  by  enforcing  mass  conservation  in  our  iterative  solution  process  so  that  the  correct  lev- 
el of  hydrostatic  pressure  is  induced  in  the  enclosed  phase.  This  is  done  in  the  following  way. 
In  order  to  satisfy  the  mass  constraint  we  require,  in  each  enclosed  phase: 


^bubble 


dx 

IS  —  cunsiani 

(2.48) 


QHy—ds  =  constant 


where  Hy  is  Jiy^  for  axisymmetric  and  y  for  planar  case,  q  is  the  density  of  the  fluid  inside 
the  enclosed  mass,  such  as  a  drop,  s  is  the  arc  length  and  s„jax  is  the  total  arc  length  of  the 
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curve  representing  the  boundary  of  the  enclosed  mass.  In  Eq.  (2.48),  a  bubble  is  used  as  the 
example;  obviously,  the  same  procedure  is  applicable  to  other  cases,  such  as  a  drop.  Since 
the  pressure  correction  formulation  only  provides  an  estimate  of  the  pressure  gradients  nec- 
essary to  satisfy  the  continuity  equation  cellwise  (i.e.,  across  control  volume  faces),  global 
conservation  information  is  not  provided  by  this  means.  Thus,  the  continuity  constraint  as 
expressed  by  the  pressure  correction  equation  establishes  the  hydrodynamic  part  of  the  pres- 
sure inside  the  enclosed  mass  consistent  with  the  flow  field,  but  provides  no  information  re- 
garding the  hydrostatic  component. 

The  dynamic  form  of  the  Young-Laplace  equation  at  the  interface  establishes  the  nor- 
mal stress  balance  as: 

A{pj  +  pc)  +  A{viscous  terms)  =  y(jt^  +  x-^  (2.49) 

where  represents  the  jump  across  the  interface  and  the  viscous  terms  take  the  appropriate 
forms  in  planar  and  axisymmetric  cases,  as  given  in  Eqs.  (2.8-2.9).  In  the  above  equation, 
subscript  d  indicates  the  hydrodynamic  and  c  the  hydrostatic  parts  of  the  pressure.  The  hy- 
drostatic pressure  must  be  such  that  it  satisfies  the  mass  constraint. 

Kang  and  Leal  (1987)  give  the  appropriate  expression  for  pc  term  as  a  time-dependent 
parameter  C(t).  This  constant  arises  when  the  pressure  at  the  interface  is  obtained  by  inte- 
grating the  component  of  the  momentum  equations  in  the  direction  tangential  to  the  interface 
in  each  phase.  This  is  easier  to  see  in  the  case  where  the  flow  on  each  side  of  the  drop  is 
obtained  separately  and  then  matched  at  the  interface  through  the  jump  condition  Eq.  (2. 12), 
as  is  the  practice  in  the  boundary-fitted  formulation  employed  by  Kang  and  Leal  (1987). 
In  related  work,  the  practice  of  rescaling  the  interface  shape  periodically  to  maintain  the  en- 
closed mass  constant  is  adopted  (Stone  and  Leal,  1989a).  In  the  work  presented  by  Nobari 
et  al.  (1996),  it  is  mentioned  that  the  mass  of  drops  is  maintained,  but  the  mechanism  by 
which  this  is  established  by  their  solver  is  not  clear.  In  recent  work  on  the  dynamics  of  menis- 
ci in  crystal  growth,  Shyy  and  Rao  (1995),  with  the  pressure-correction  methodology  used 
here,  maintain  the  volume  of  the  meniscus  by  employing  an  iterative  procedure  to  obtain  the 
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proper  value  of  the  hydrostatic  pressure,  consistent  with  Eq.  (2.49)  as  well  as  mass  conserva- 
tion. However,  the  computations  use  the  boundary-fitted  technique  and  the  implementation 
of  the  correction  of  mass  conservation  need  to  be  modified  in  the  present  context. 

A  method  for  correcting  the  hydrostatic  pressure  in  the  enclosed  mass  needs  to  be  de- 
vised so  that  it  does  not  lose  or  gain  mass.  To  achieve  this  goal,  we  devise  an  indirect  method 
of  establishing  mass  conservation  similar  to  that  adopted  by  Kang  and  Leal  ( 1 987).  The  de- 
sired result  in  enforcing  mass  conservation  is  to  establish  an  additive  constant  to  the  pressure 
inside  the  enclosed  mass.  The  stress  balance  equation  in  Eq.  (2.49)  then  becomes: 

A(p^  +  Pc  +  dpc)  +  A{viscous  terms)  =  yix^  +      +  dx)  (2.50) 

The  correction  to  the  pressure,  dpc,  is  responsible  for  changes  in  the  drop  shape  repre- 
sented by  dx  so  that  upon  integrating  as  in  Eq.  (2.49)  over  this  new  shape,  the  volume  is  con- 
served. In  our  procedure,  during  the  iterative  process  of  establishing  the  flow  field,  we  im- 
pose the  mass  constraint  on  the  interface  motion  by  establishing  the  shape  correction  dx  and 
inducing  the  pressure  correct  increment  dpc.  At  each  iteration,  this  shape  correction  is  per- 
formed after  the  interface  has  been  advected  to  its  new  position  based  on  the  kinematic  condi- 
tion. Thus  the  fresh  iterate  of  the  j:-coordinate  of  interfacial  marker  position  is  obtained  by 

(X^'''''^^)p  =  X/  +  dtVnn, 

)c  -  (^f       )p  +  o„nx  (2.51) 

where  the  subscripts  p  and  c  represent  the  predicted  and  corrected  values,  the  subscript  £  is 
the  interfacial  marker  index,  the  superscripts  and  k  represent  the  iteration  number  and  time 
step,  respectively.  Similar  expressions  are  obtained  for  the  ^'-coordinate. 

The  shape  correction  in  the  normal  direction  6^  is  performed  in  an  iterative  fashion  by 
means  of  a  bisection  method,  so  that  the  mass  of  the  drop  as  computed  in  Eq.  (2.48)  remains 
a  constant.  The  value  of  d„  is  obtained  from  the  knowledge  of  the  arc  length  measure  of  the 
interface  and  the  estimate  of  the  mass  loss.  The  bisection  method  used  to  obtain  the  value 
of  (5„  converges  in  a  few  iteration  (typically  <  5).  The  entire  mass  conservation  measure  is 
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enclosed  in  the  outer  iteration  loop  for  the  entire  solution  procedure,  including  interface 
movement,  flow  field  computation  and  mass  constraint  imposition.  The  procedure  is  then 
required  to  achieve  a  preset  tolerance  for  convergence  levels  of  the  flow  field,  interface  shape 
and  enclosed  mass.  It  should  be  emphasized  that  the  present  procedure  is  of  a  predictor-cor- 
rector type,  and,  at  convergence  should  yield  a  drop  shape  independent  of  the  manner  dpc 
is  imposed. 

1.6.  Summary 

In  this  chapter  we  have  described  the  fundamental  elements  of  the  computational  meth- 
odology called  ELAFINT.  The  numerical  framework  developed  here  have  several  improve- 
ments over  the  previous  version  of  the  code  that  can  be  summarized  as  follows: 

( 1 )  .  A  collocated  variable  approach  is  adopted.  The  stability  and  accuracy  of  the  algorithm 

in  computing  flows  of  comparable  complexity  to  those  handled  by  the  staggered  grid 
approach  was  established.  The  collocated  layout  considerably  simplifies  the  algorithm. 

(2)  .  The  cut  cell  technique  is  facilitated  by  the  collocated  grid  layout,  and  is  rendered  less 

tedious.  While  previously  the  cut  cell  technique  had  to  contend  with  three  types  of  con- 
trol volumes,  in  the  present  case  only  one  control  volume  is  used. 

(3)  .  The  fluid-fluid  interface  is  modeled  using  the  inmiersed  boundary  technique.  The  cut- 

cell  approach  was  used  in  previous  work  on  solidification  for  the  moving  boundary  cal- 
culations. For  fluid-fluid  interfaces  however,  this  method  proves  to  be  sensitive  to 
noise  and  difficult  to  stabilize.  This  is  because  the  cut-cell  method  employs  no  dissipa- 
tion, so  that  the  boundary  condition  application  is  exact.  In  other  moving  boundary 
treatments  based  on  a  fixed  grid,  some  smearing  is  added  near  the  discontinuity.  This 
is  the  case  in  the  immersed  boundary  technique  (Unverdi  and  Tryggvason,  1992),  the 
VOF  method  (Kothe  et  al.,  1992)  and  the  level-set  method  (Sussman  et  al.,  1992).  The 
choice  of  the  immersed  boundary  technique  is  made  due  to  its  compatibility  with  ex- 
plicit interface  tracking  where  ambiguity  in  regard  to  interface  location  is  avoid. 


CHAPTER  3 

ASSESSMENT  OF  COMPUTATIONAL  TECHNIQUE 

In  the  preceding  chapter,  the  elements  of  the  numerical  framework  developed  for  han- 
dling moving  boundary  problems  involving  fluid-solid  interaction  as  well  as  multi-fluid 
flows  are  described.  In  this  chapter,  the  focus  is  to  illustrate  the  accuracy  and  versatility  af- 
forded by  the  combination  of  the  cut-cell  method  and  the  immersed  boundary  technique. 
Towards  this  objective,  we  present  several  designated  cases  to  evaluate  the  performance  of 
the  cut-cell  method  in  the  current  coUocated-variable  framework.  This  cut-cell  approach 
was  tested  in  previous  work  (Shyy  et  al.,  1996;  Udaykumar  et  al.,  1996)  for  a  staggered  vari- 
able arrangement.  Results  given  by  the  immersed  boundary  technique  in  the  present  finite- 
volume,  pressure-based  framework  are  compared  to  known  solutions  of  viscous  drop  dy- 
namics. For  this  particular  application,  reliable  experimental  and  numerical  results  exist,  at 
least  for  highly  viscous  flow,  with  which  we  compare  our  results.  Finally,  we  put  the  cut-cell 
and  immersed  boundary  representations  together  to  enable  the  simulation  of  drop  dynamics 
in  a  constricted  tube.  This  last  aspect  demonstrates  the  ability  of  the  present  method  to  handle 
multi-fluid  flows  in  arbitrary  geometries  under  viscosity-  or  inertia-dominated  flow  condi- 
tions. 

3.1.  Flows  With  Fixed  Arbitrary-Shaped  Internal  Boundaries 

One  objective  of  this  work  is  to  develop  a  facility  for  computing  complex  geometry 
flows  on  Cartesian  grids.  This  feature  will  remove  restrictions  regarding  the  flow  domains 
in  which  the  moving  boundary  simulations  can  be  performed,  without  having  to  resort  to 
body-fitted  grids  or  unstructured  grid  layouts. 

We  first  compare  solutions  for  a  stationary  interface  with  those  from  well-tested  solu- 
tion procedures  using  body-fitted  coordinates  (Shyy,  1 994).  This  is  one  of  several  exercises 
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(a) 


Figure  3.1.  Comparison  of  a  deformed  cavity  at  Re=  1000,  for  results  using  the 
cut-cell  technique  on  the  collocated  grid  with  the  solution  on  a  boundary-fitted 
staggered  grid,  (a)  Streamline  contour  for  the  collocated  grid,  (b)  u-velocity 
along  vertical  centerline  for  the  present  collocated  grid  case  (solid  line)  and  the 
staggered  grid  case  (open  circles),  (c)  v- velocity  comparison  along  vertical 
centerline. 
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performed  to  validate  the  conservation  and  consistency  characteristics  of  the  cut-cell  meth- 
od used  to  treat  the  complex  geometry.  The  interface  is  sufficiently  deformed  that  all  the 
possible  types  of  cut  cells  are  encountered.  The  square  cavity  is  a  frequently  adopted  test  i 

bed  for  numerical  experiments  for  incompressible  flows  due  to  available  benchmarked  solu-  | 

■J 

tions  for  comparisons.  We  deform  the  base  of  the  cavity,  the  amplitude  being  10%  of  the 
base.  A  121x121  Cartesian  grid  is  employed.  We  first  present  the  results  for  a  driven  cavity 
flow,  where  the  top  wall  of  the  cavity  is  pulled  at  velocity  U=l  corresponding  to  a  Reynolds 
number  of  1 000.  The  results  from  the  present  method  are  compared  with  a  previously  bench- 
marked  body-fitted  staggered  grid  formulation  (Thakur  et  al . ,  1 996)  with  the  same  grid  size . 
The  results  are  shown  in  Figure  3.1.  The  streamline  pattern  is  shown  in  Figure  3.1(a).  A 
quantitative  comparison  can  be  seen  in  Figures  3 . 1  (b)  and  (c)  where  the  centerline  velocities 
are  plotted.  As  seen  in  this  figure,  the  centerline  profiles  obtained  from  the  cut-cell  tech- 
nique using  the  collocated  grid  (solid  line)  agree  well  with  the  body-fitted  staggered  grid 
calculations  (open  circles). 

We  next  demonstrate  the  ability  of  the  complex  geometry  treatment  to  simulate  more 

-I 

complex  flow  situations.  In  Figure  3.2,  the  streamlines  and  pressure  contours  for  a  driven 
cavity  with  top  and  bottom  moving  walls,  and  enclosing  a  semicircular  cavity  at  the  bottom 
wall  are  shown.  The  smaller  cavity  is  indicated  with  dotted  lines  in  Figure  3.2,  which  repre- 
sents a  solid  boundary  with  zero  internal  thickness.  This  configuration  gives  an  indication 
of  the  versatility  of  the  present  algorithm. 

Figures  3.2(a)  and  (b)  show  the  streamlines  and  the  pressure  contours  in  the  two  cavi- 
ties, respectively.  As  can  be  seen  from  the  figure,  for  the  Reynolds  number  of  1000,  the  mul- 
tiple recirculation  zones,  including  the  one  in  the  smaller  cavity  are  well  captured.  It  would 
be  a  significant  task  to  fit  a  good  boundary-fitted  grid  to  this  configuration.  An  added  advan- 
tage that  may  be  realized  is  that  the  grid  can  be  refined  easily  in  the  present  Cartesian  frame- 
work without  consideration  of  the  solid  boundary  shape  or  the  number  of  solid  boundaries 
present. 
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Figure  3.2.  Double-Driven  cavity  flow  with  obstruction.  Re=1000,  121x121 
grid.  Dotted  line  shows  the  shape  of  the  enclosed  cavity,  (a)  Contour  plot  for 
streamline,  (b)  For  pressure. 
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Figure  3.3  shows  the  results  of  a  natural  convection  flow  in  a  square  cavity  with  multi- 
ple objects  in  the  flow  field.  The  results  are  for  a  flow  with  a  Prandtl  number  of  1  and  Ray- 
leigh  number  of  10^.  The  sidewalls  of  the  cavity  are  heated,  the  non-dimensionalized  tem- 
perature of  the  left  wall  maintained  at  0  and  that  of  the  right  wall  at  1.  The  top  and  bottom 
walls  are  adiabatic.  Each  of  the  cylinders  is  maintained  at  a  constant  temperature  of  0. 

In  Figure  3.3(a)  we  show  the  isotherms  and  in  Figure  3.3(b)  the  velocity  vectors.  The 
skew  symmetry  of  the  temperature  field,  obtained  in  such  a  cavity  without  obstacles,  is  bro- 
ken by  the  placement  of  the  five  shapes.  In  particular,  the  temperature  gradients  in  the  right 
bottom  comer  appear  to  be  enhanced  by  the  presence  of  the  cylinder  in  that  location  while 
those  in  the  lower  left  comer  are  diminished.  This  would  obviously  have  a  significant  impact 
on  the  heat  transfer  within  the  cavity  shown.  This  type  of  a  facility  would  be  useful  in  com- 
puting heat  transfer  around  complex  geometries,  and  can  be  combined  with  phase  change 
as  demonstrated  in  Udaykumar  et  al.  (1996). 

No  restriction  exists  on  the  number  or  shape  of  the  obstructions  that  can  be  inserted  in 
the  domain.  Since  all  the  measures  dealing  with  the  interior  solid  boundaries  are  handled 
by  means  of  one-dimensional  arrays,  the  expense  of  incorporating  intemal  boundaries  varies 
linearly  with  the  number  of  objects.  No  grid  generation  measures  are  required  when  a  new 
object  is  introduced,  thus  objects  can  be  added  or  deleted  without  affecting  the  existing  ge- 
ometry globally. 

Finally,  a  transient  calculation  is  performed  in  an  open  channel  configuration.  The  case 
shown  in  Figure  3.4  demonstrates  the  usefulness  of  this  procedure  for  the  time-dependent 
computation.  Here,  the  inlet  is  at  the  left  of  the  domain.  The  Reynolds  number  in  this  case 
is  100.  A  121x81  grid  is  used.  The  development  of  the  recirculation  zone  behind  the  cylinder 
is  shown  until  a  steady-state  is  reached  at  a  non-dimensional  time  t=3.5. 

In  Figure  3.4(e),  we  show  the  streamlines  obtained  for  this  geometry  from  a  steady- 
state  calculation  on  a  well-tested  boundary-fitted  staggered  grid  (Thakur  et  al.,  1996).  The 
flow  fields  obtained  from  the  two  methods  are  identical.  In  particular,  the  recirculation  zone 


Figure3.3.  Natural  convection  flow  in  a  cavity.  Pr=l,Ra=10^,  121x121  grid, 
(a)  Contour  plot  for  temperature,  (b)  Velocity  vectors. 
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(a)  t*=0.5 


Figure  3.4.  Open  channel  flow,  transient  calculation.  Re=100,  121x81  grid, 
(a)  to  (d)  Streamline  plot  showing  development  of  the  steady  recirculating 
flow  behind  the  cylindrical  obstruction  computed  using  the  cut-cell  method, 
(e)  Steady-state  streamlines  from  the  body-fitted  calculation. 

in  each  case  extends  up  to  x=3.61  in  the  rear  of  the  cylindrical  obstruction.  The  slight  differ- 
ences in  the  streamlines  between  (d)  and  (e)  are  due  to  the  differences  in  data  representation 
on  the  Cartesian  and  curvilinear  grids  which  is  interpreted  differently  by  the  contouring  fa- 
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cility.  It  is  also  noted  that  the  cylindrical  obstruction  has  been  shown  in  this  figure  by  means 
of  a  phase  fraction  contour.  The  step-like  nature  of  the  obstacle  is  a  limitation  of  the  contour- 
ing procedure.  As  far  as  the  solver  itself  is  concerned,  the  surface  is  represented  smoothly 
by  means  of  the  cut-cell  technique. 

3.2.  Moving  Boundary  Flow  Problems 

3.2. 1 .  Deformation  and  Recovery  of  Viscous  Drops  in  Extensional  Flows 

We  now  present  results  for  the  drop  deformation  and  recovery  studies.  Consider  the 
motion  of  an  incompressible,  neutrally  buoyant  drop  in  a  straining  flow  described  by  the 
Stokes  equations.  The  non-dimensionalization  of  the  governing  equations  is  performed  us- 
ing the  reference  quantities  as  described  below. 

We  define  the  viscosity  ratio  between  drop  and  suspending  fluids  as  A  =  Hi/Ho  ■  The  char- 
acteristic velocity,  pressure,  viscosity  and  length  are  the  following:     =  y/fic<  Pc  =  (f^c 
€c),/^c  =  Mo  is  the  viscosity  of  the  suspending  fluid,  (,^  =  a  =  the  radius  of  undeformed  circular 
shape  and  y  is  the  surface  tension.  Then,  the  non-dimensionalized  governing  equations  for 
creeping  flows  become 
Continuity  equation: 

V  •  V  =  0  (3.1) 
Momentum  equation  for  suspending  fluid: 

0  =  -  Vp  +  V^V  +     xnd{x  -  (Xf))ds 

J  (3.2) 

for  the  drop: 

0  =  -Vp  +  AV^V  +     xnd(x  -  {x^))ds 

J  (3.3) 

In  the  above  equations,  V  and  p  represent  the  velocity  and  pressure,  respectively,  n  is 
the  unit  normal  vector,  defined  as  positive  when  pointing  outward  from  the  drop  phase  to 
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the  suspending  fluid  phase,  x  =  V  •  n  represents  the  local  mean  curvature  of  the  interface, 
as  defined  in  Chapter  2.2.2;  6  denotes  the  delta  function,  3c  and  (3c)^  are  the  location  of  points 
in  the  flow  field  and  on  the  interface.  The  integration  term  shown  in  Eqs.  (3.2)  and  (3.3) 
represents  the  Eulerian  contributions  of  the  normal  stress  discontinuity  due  to  the  capillary 
force  evaluated  on  the  interface.  This  force  is  transmitted  from  the  interface  to  the  fluid  by 
use  of  the  Immersed  Boundary  Technique  originated  by  Peskin  (1977).  The  detailed  descrip- 
tion of  the  numerical  implementation  of  the  force  term  is  given  in  Chapter  2.4.2. 
The  far  field  condition  for  planar  flow  is 
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and  for  axisymmetric  flow. 
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where  the  capillary  number  Ca  which  represents  the  ratio  of  viscous  force  to  surface  tension 
is  given  by  (f^g  Ga/y).  G  is  the  strain  rate  of  the  imposed  flow.  The  discretization  of  the  above 
equations  is  described  in  Chapter  2.2.4. 

The  flow  condition  and  grid  layout  is  indicated  in  Figure  3.5.  The  grid  distribution  is 
nonuniform;  in  the  far  field,  i.e.,  away  from  the  inner  region  shown,  the  grids  are  stretched 
linearly  toward  the  boundaries  of  the  domain,  while  in  the  inner  region  a  uniform  fine  grid 
is  used.  This  yields  better  resolution  in  the  desired  region  of  the  flow  field,  i.e.,  in  the  vicinity 
of  the  interface.  The  flow  field  is  imposed  as  shown,  thereby  exploiting  the  symmetry  of 
the  problem.  The  incoming  flow  at  the  top  and  outgoing  flow  at  the  right  are  held  constant 
as  specified  by  Eq.  (3.4)  and  Eq.  (3.5)  for  planar  and  axisymmetric  cases,  respectively.  The 
distribution  of  markers  along  the  interface  is  maintained  at  a  constant  arc  length  spacing. 
The  arc  length  spacing  used  in  all  computations  presented  here  is  0.8h,  which  was  found  by 
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Regio  1  of  coarse  grid 


Axis  of  symmetry 


Region  of  fine  grid 


Outgoing  flow 


Axis  of  symmetry 


Location  of  droplet 


Figure  3.5.  Illustration  of  domain  of  computation  and  boundary  conditions 
for  drop  deformation  simulation. 


experimentation  to  provide  adequate  resolution  and  accuracy.  The  redistribution  along  the 
interfacial  curve  is  performed  at  every  time  step  to  prevent  undesirable  depletion  or  accu- 
mulation of  markers.  The  dimensions  L  and  W  are,  respectively,  defined  as  the  deformed 
extents  of  the  drop  in  the  x-  and  y-  directions. 

At  the  outset,  a  suitable  grid  is  determined  for  the  moving  boundary  simulations.  A  grid 
independence  study  was  first  performed  to  obtain  the  steady-state  shapes  of  the  drop  in  ex- 
tensional  flow  at  a  capillary  number  of  0. 1 ,  which  is  slightly  below  the  critical  capillary  num- 
ber. In  Figure  3.6,  we  compare  the  results  obtained  from  calculations  on  an  81x81  and 
161x161  grid  with  and  without  enforcing  mass  conservation  condition.  As  can  be  seen,  the 
results  from  the  two  grids  with  mass  conservation  imposed  are  in  close  agreement.  We 
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Figure  3.6.  Grid  independence  test  for  Ca=0. 1 ,  X,= 1 .  The  final  steady-state  siiapes 

are  shown  for  (a)        :  81x81  grid  without  mass  conservation  (b)  —  :  161x161 

grid  without  mass  conservation  (c)  :  81x81  grid  with  mass  conservation  (d) 

-.-  :  161x161  grid  with  mass  conservation. 


henceforth  retain  the  same  grid  density  as  that  corresponding  to  the  81x81  grid  in  the  present 
case.  The  number  of  grid  points  will  therefore  scale  with  the  domain  size. 

The  moving  boundary  capability  using  the  immersed  boundary  technique  is  first  tested 
against  the  results  of  Stone  and  Leal  (1989a,  b)  for  steady  axisymmetric  drop  shapes  under 
an  uniaxial  extensional  flow.  For  a  range  of  capillary  numbers  below  the  critical  capillary 
number  Cctr  the  drop  deformation  culminates  in  a  steady-state  shape.  Stone  and  Leal 
(1989a)  have  proposed  that  the  deformation  parameter  D=(L-W)/(L+W)  is  an  indicator  of 
the  final  bubble  shape  at  steady-state,  where  L  is  the  major  radius  and  W  the  minor  radius 
of  the  ellipsoidal  drop  as  shown  in  Figure  3.5.  This  quantity  therefore  provides  an  estimate 
of  the  aspect  ratio  of  the  steady-state  shape. 
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In  Figure  3.7,  we  compare  the  D  vs  Ca  curves  obtained  from  our  calculations  with  those 
of  Stone  and  Leal  (1989a,b)  for  viscosity  ratio  A =7.  The  grid  in  this  case  was  61x61  in  the 
inner  region  of  dimensions  3x3  units.  The  arc  length  spacing  of  the  interfacial  points,  follow- 
ing testing  to  determine  accuracy,  was  fixed  to  be  0.8h  in  this  and  following  calculations. 
The  number  of  points  on  the  interface  was  45  at  steady  state.  As  can  be  seen  from  the  figure, 
our  results  are  in  good  agreement  with  those  of  Stone  and  Leal  (1989a).  Furthermore,  the 
critical  capillary  number  of  Cacr=0.12  obtained  in  the  present  case  also  corresponds  to  the 
value  in  Stone  and  Leal  (1989a). 

Figures  3.7(b)  and  (c)  show  the  flow  field  and  pressure  contours  in  the  case  of  Ca=0.1. 
It  is  clear  from  the  velocity  vector  plot  that  the  flow  at  steady-state  is  tangent  to  the  surface 
of  the  drop  and  the  single-domain  treatment  of  the  interfacial  forces  yields  the  desired  shape 
and  flow  field  by  the  balance  of  normal  stresses  in  Eq.  (2. 12).  The  large  pressure  gradient 
across  the  interface  is  also  seen  clearly  in  this  figure. 

We  next  study  the  unsteady  deformation  and  recovery  characteristics  of  drops  under  im- 
posed straining  flow  with  Ca  >  Cocr-  The  inner  region  of  fine  grid  in  this  case  is  6x6  units. 
The  domain  of  computation  is  8x8  units  in  size.  In  the  inner  region  121x121  grids  are  used, 
while  in  the  coarse  grid  region  20  points  are  used  along  each  direction  and  the  grid  is  linearly 
stretched.  In  Figure  3.8(a),  we  show  the  shapes  of  the  drop  at  different  times  for  the  Capillary 
number  of  0. 14.  In  this  case,  the  drop  deforms  to  a  maximum  extent  of  5.4  times  the  initial 
radius  at  non-dimensional  time  t=35.  The  number  of  points  on  the  interface  increases  from 
27  initially  to  140  at  the  final  stage. 

To  study  the  recovery  of  the  drop,  we  turn  off  the  flow  at  the  stages  corresponding  to 
the  shapes  at  t=28  and  t=35  and  track  the  recovery  path  in  the  two  cases.  The  results  are 
presented  in  Figures  3.8(b)  and  (c).  As  can  be  seen,  for  the  smaller  deformation  case  the 
interface  recovers  back  to  the  spherical  drop  shape,  while  for  the  larger  deformation  case, 
the  recovery  process  results  in  the  impending  breakup  stage  shown.  Although  the  further 
development  of  the  drop  after  this  stage  can  be  handled  by  the  code,  we  stop  the  computations 
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Figure  3.7.  (a)  Comparison  of  D=(L-W)/(L+W)  vs.  Ca  obtained  from  the  pres- 
ent method  (solid  line)  with  the  results  (dashed  line)  of  Stone  and  Leal  (1989a). 
The  critical  capillary  number  is  approximately  0. 12  in  both  studies  for  l.=l.  (b) 
Velocity  vectors  at  steady-state  for  Ca=0. 1 .  (c)  Pressure  contours  at  steady-state 
for  Ca=0.1. 


53 


(a) 


t  =  0,  7,21,28,  35 


(b) 


t  =  0,5,  10,  25 


(c) 


t  =  0,  10,  27,31 


Figure  3.8.  Deformation  and  relaxation  of  the  drop  for  capillary  number  of  0.14 
(critical  capillary  number  is  0. 12).  (a)  Deformation  stage,  shapes  shown  at  t=0,  7, 
21 , 28  and  35.  At  t=35,  the  deformation  is  5.4  times  the  original  undeformed  radius, 
(b)  Recovery  of  the  drop  from  shape  shown  at  t=28  in  (a),  (c)  Recovery  of  the  drop 
from  shape  shown  at  t=35  in  (a). 


at  this  stage.  From  Figure  3.8(b),  it  is  seen  that  the  drop  returns  to  a  sphere  of  radius  1 ,  which 
demonstrates  that  the  mass  conservation  procedure  implemented  holds  well  in  the  course  of 
the  large  deformations  experienced  by  the  drop.  It  may  be  surmised  from  Figure  3.8(c)  that 
a  small  satellite  drop  is  formed  along  with  two  larger  drops  in  the  recovery  from  the  large 
deformation.  This  behavior  of  the  drops  with  extent  of  deformation  for  the  same  capillary 
number  is  in  qualitative  agreement  with  the  results  of  Stone  and  Leal  (1989b). 
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In  Figure  3.9,  we  show  the  flow  fields  resulting  from  the  computations.  In  Figure 
3.9(a),  the  velocity  vector  is  at  the  top  and  the  pressure  contour  is  at  the  bottom  for  a  drop 
deformed  at  t=35  and  Ca=0. 14.  As  can  be  seen,  the  flow  in  this  case  tends  to  draw  the  drop 
out  into  a  filament.  Figure  3.9(b)  shows  the  velocity  vectors  for  the  recovered  shape  after 
release  from  a  small  deformation.  The  velocity  magnitudes  are  very  low  (O(IO^))  in  this 
nearly  quiescent  state  of  the  flow  field.  The  corresponding  pressure  contours  are  shown  un- 
derneath. It  can  be  seen  that  the  pressure  levels  inside  the  drop  are  uniform  and  equal  to 
approximately  1.96  in  value,  which  is  close  to  the  expected  static  value  for  the  given  initial 
drop  shape. 

In  Figure  3.9(c),  the  velocity  vectors  and  corresponding  pressure  contours  inside  a  drop 
corresponding  to  the  recovery  case  shown  at  t=31  in  Figure  3.8(c)  are  shown.  The  impend- 
ing breakup  of  the  drop  is  evident  from  the  large  velocity  and  pressure  values  observed  in 
this  figure.  The  capillary  induced  flow  which  results  in  the  breakup  of  the  drop  is  clearly 
seen  in  regions  of  large  negative  curvature  in  Figure  3.9(c).  Due  to  the  lack  of  resolution 
in  the  region  of  breakup,  the  values  in  this  figure  are  not  of  significance.  It  is  noteworthy 
that  the  collocated  grid  solution  framework  holds  up  in  the  presence  of  the  large  pressure 
gradient  seen  in  this  case  (maximum  pressure  value=24,  minimum=-0.8). 

Figure  3. 10  presents  the  recovery  characteristics  of  viscous  drops  in  uniaxial  extension- 
al  flow  for  the  viscosity  ratios  A=i  and  10.  In  Figure  3.10(a),  we  plot  the  deformation  of 
a  drop  for  the  case  of  Ca=0.13,  which  is  slightly  above  the  critical  capillary  number,  and 
X=l,  i.e.,  the  viscosities  of  the  drop  and  the  suspending  fluid  are  the  same.  The  drop  is  de- 
formed to  a  certain  L  and  a  step  change  in  the  flow  to  Ca=0.065  is  made.  The  step  change 
in  the  imposed  flow  strength  is  made  at  three  different  points  along  the  deformation  curve. 
The  shapes  along  the  deformation  and  recovery  curves  are  also  shown  in  the  inset  boxes. 
These  curves  seek  to  demonstrate  the  accuracy  of  the  present  calculations  for  long  periods 
and  large  magnitudes  of  drop  deformation.  The  results  here  agree  very  well  with  those  of 
Stone  and  Leal  (1989b).  In  particular,  the  fate  of  the  recovering  drops  from  three  different 
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Figure  3.9.  Velocity  vectors  and  corresponding  pressure  contours  plots,  (a)  deformed 
drop  at  the  time  t=35  in  Figure  3.8(a);  (b)  drop  at  a  relaxation  time  t=25  in  Figure  3.8(b); 
(c)  drop  at  a  relaxation  time  t=31  in  Figure  3.8(c). 
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initial  deformations,  following  a  step  change  in  the  flow  strength  match  with  the  predictions 
in  Stone  and  Leal  (1989b). 

The  initial  large  deformation  of  the  drop  due  to  stretching  is  evident  in  the  inset  marked 
I.  When  the  released  drop  takes  path  II,  i.e.,  for  a  L  value  before  release  of  2. 79,  the  drop 
recovers  to  a  steady-state  shape  corresponding  to  Ca=0.065.  For  a  value  of  L=3.07,  no  re- 
covery is  possible.  In  this  case,  as  can  be  seen  in  the  recovery  curve,  the  drop  initially  shrinks. 
However,  during  this  process  a  waist  develops  at  the  equator  of  the  drop  and  results  in  the 
extension  of  the  drop  as  the  flow  induced  by  the  capillary  forces  accompanying  this  negative 
curvature  drives  the  fluid  outward,  as  seen  in  the  extreme  case  in  Figure  3.9(c).  This  situation 
is  unstable  and  leads  to  the  further  development  of  the  waist  and  subsequent  breakup  of  the 
drop.  Thus,  drop  breakup  in  such  flows  is  dependent  on  the  competition  between  the  im- 
posed flow  fields  and  the  flow  induced  by  capillarity.  Therefore,  although  a  waist  exists  in 
the  deformation  curve  marked  I,  in  this  situation  the  flows  due  to  stretching  and  capillarity 
reinforce  each  other,  while  in  the  recovery  phase  they  oppose  each  other. 

The  result  of  the  recovery  process  depends  on  the  strength  of  the  destabilizing  flow  in- 
duced by  the  capillary  forces  acting  at  the  waist,  which  depends  on  the  curvature  at  the  waist. 
In  Figure  3.10(b)  for  example,  it  is  seen  that  the  drop  in  case  IV,  which  represents  L=3.24 
before  release,  rapidly  heads  towards  breakup.  This  is  in  contrast  to  the  other  cases,  where 
there  is  an  initial  tendency  to  respond  to  the  imposed  flow  field  before  capillarity  takes  over. 
It  is  therefore  clear  that  questions  regarding  recovery,  number  of  drops  resulting  from  break- 
up, sizes  of  drops  and  so  on  can  only  be  answered  in  the  context  of  the  flow  fields  imposed, 
relative  properties  of  drops  and  suspending  fluids  and  extent  of  deformation. 

The  rates  of  deformation  obtained  from  our  calculations  agree  very  well  with  those  in 
Stone  and  Leal  (1989b).  Some  differences  in  times  of  deformation  exist  between  the  two 
calculations  in  the  initial  phase  of  the  deformation.  These  may  be  due  to  the  differences  in 
initial  conditions  or  capillary  numbers  chosen  for  the  deformation.  A  notable  point  is  that 
the  time  scale  of  deformation  of  drops  in  viscous  extensional  flows  is  expected  to  be  propor- 
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Figure  3.10.  Recovery  curves  for  viscous  drop  subjected  to  uniaxial  extensional 
flow  as  a  test  of  the  accuracy  of  the  method  for  large  deformations  and  property 
jumps,  (a)  The  recovery  curve  for  X=l.  (b)  Recovery  curve  for  X=  10. 
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tional  to  (7+A)  (Stone  and  Leal  1989a).  A  comparison  between  Figures  3.10(a)  and  (b)  re- 
veals the  validity  of  this  scaling. 

3.2.2.  Deformation  of  Drops  in  Extensional  Flows  with  Inertia  Effects 

In  order  to  demonstrate  the  ability  of  the  algorithm  to  deal  with  convective  drop  dynam- 
ics, we  investigate  the  effects  of  an  imposed  uniaxial  extensional  flow.  The  scales  chosen 
for  this  problem  are  the  same  as  the  ones  for  the  viscous  drop  presented  in  the  previous  sec- 
tion. However,  in  this  case,  with  the  inclusion  of  the  convection  terms,  the  non-dimensional- 
ized  governing  equations  are 
Continuity  equation: 

V  •  V  =  0  (3.6) 
Momentum  equation  for  suspending  fluid: 


+  V(VV)  =  -Vp  +  -^V^V  +  4-  I  xnd(x  -  ix^))ds 


(3.7) 


for  the  drop: 


9X  +  V(Vy)  =  -Vp  +  j-V^V  +  -ji-  I  xndix  -  (Xf))ds 

^  (3.8) 


1 


Under  convective  conditions  the  relevant  non-dimensional  parameters  that  arise  are  the 
Reynolds  number  {Re=2Qc(G£cHc/Mc^  based  on  the  drop  radius)  and  Weber  number 
(We=2Qc(G€c)'^icA^)-  The  former  provides  a  measure  of  the  relative  dominance  of  inertia 
to  viscous  effects,  while  the  latter  is  the  ratio  of  inertia  to  capillary  forces.  In  the  results  pres- 
ented the  imposed  flow  corresponds  to  a  Reynolds  number  of  10.  We  attempt  to  obtain  the 
critical  Weber  number  for  a  drop  under  such  an  imposed  linear  flow  condition.  Kang  and 
Leal  (1987)  have  obtained  the  critical  Weber  number  of  a  bubble  at  this  Reynolds  number 
as  lying  between  0.9  and  1 .0.  In  that  work,  the  density  and  viscosity  inside  the  bubble  was 
considered  to  be  negligible  so  that  the  pressure  inside  the  bubble  remains  a  constant,  and  the 
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interior  of  the  bubble  is  passive,  in  that  it  has  no  effect  on  the  deformation  characteristics. 
Under  this  assumption  they  fit  a  boundary-conforming  grid  to  the  exterior  of  the  bubble  and 
obtained  the  transient  deformation  characteristics  of  the  bubble  surface.  This  treatment  is 
a  limitation  imposed  by  the  boundary-fitted  grid  approach.  If  calculations  were  to  be  per- 
formed inside  the  bubble  also,  grid-generation  would  become  complicated.  In  our  frame- 
work this  limitation  does  not  exist.  On  the  other  hand,  it  is  not  a  simple  matter  to  render  the 
flow  field  inside  the  bubble  entirely  passive,  since  the  present  method  is  a  single-domain 
method.  This  rendering  of  a  part  of  the  domain  passive  was  accomplished  in  previous  cal- 
culations involving  stadonary  and  moving  solid  boundaries,  where  there  is  no  flow  in  the 
solid  region.  In  the  case  of  two-fluid  interfaces  this  is  not  as  straightforward  to  do  due  to 
the  fact  that  the  interface  velocity  is  governed  by  the  kinematic  condition  at  the  boundary, 
and  the  boundary  velocity  is  obtained  from  a  single-domain  formulation,  Eq.  (2.44).  Thus 
information  regarding  the  velocity  field  is  required  from  both  sides  of  the  interface,  which 
means  that  the  flow  field  needs  to  be  computed  in  both  phases.  Therefore,  in  performing  our 
calculations  in  this  secfion  we  will  only  be  able  to  obtain  a  qualitative  comparison  with  Kang 
and  Leal  (1987),  and  we  will  support  our  result  with  expectations  based  on  the  physics  of 
the  problem.  The  flow  configuration  remains  the  same  as  shown  in  Figure  3.5.  The  density 
and  viscosity  ratios  used  are  ^  =  Qi/Qq  =  0.02  and  X  =  ^li/^lo  =  0.02,  respectively. 

In  Kang  and  Leal  (1987),  it  was  shown  that  the  deformation  of  the  drop,  in  particular 
the  availability  of  a  steady-state  solution  for  a  given  Weber  number  was  dependent  on  the 
initial  state  of  the  drop  in  the  flow  field.  Although  for  certain  Weber  numbers  close  to  but 
below  critical,  gradually  increasing  the  flow  strength  from  a  lower  Weber  number  yielded 
steady-state  solutions,  imposing  a  step  change  in  the  flow  field  could  drive  the  drop  toward 
unsteady,  continuous  stretching.  Such  behavior  was  noticed  in  our  simulations  also.  There- 
fore, for  the  cases  presented  here,  we  obtained  the  flow  at  a  given  We,  under  a  particular 
combination  of  density  and  viscosity,  starting  from  the  flow  field  corresponding  to  a  value 
of  We  0. 1  lower.  The  curves  in  Figure  3. 1 1(a)  show  the  variafion  of  D  with  We  for  the  im- 
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posed  flow  at  Re=10.  To  start  this  curves,  we  first  obtained  the  flow  field  at  We=0.5  by  im- 
posing a  flow  field  instantaneously.  This  led  to  an  oscillation  of  the  drop  due  to  the  compeU- 
tion  between  inertia  and  capillarity  eff'ects.  This  type  of  behavior  is  also  reported  in 
simulations  by  Mashayek  and  Ashgriz  (1996).  For  example,  ^=0.02  and  A=0.02,  for  upon 
imposition  of  the  flow  field  the  drop  at  We=0.5  stretches  rapidly  from  its  initial  circular  shape 
to  a  value  along  the  x-axis  of  L=1.4,  and  subsequently  relaxes  in  an  oscillatory  manner  to 
a  steady-state  value  of  L= 1.17.  Stepping  up  the  Weber  numbers  in  increments  of  0. 1  thereaf- 
ter yields  the  curve  shown  in  Figure  3. 1 1  (a).  For  ^=0.02  and  A =0.02,  and  at  a  Weber  number 
of  1.2,  a  steady-state  did  not  exist.  If  these  fluid  properties  are  changed,  then  the  critical 
Weber  number  changes  accordingly. 

A  closer  examination  of  the  case  of  ^=0.02  and  X=0.02  can  be  seen  from  the  curve  of 
Lvst  shown  in  Figure  3. 1 1(b).  There  the  deformation  histories  for  We=0. 7, 1.1  and  1.2  are 
plotted.  As  can  be  seen,  these  curves  show  the  oscillatory  relaxations  mentioned  above. 
However,  in  the  case  of  We=1.2,  the  oscillations  carry  the  interface  away  from  a  steady-state 
and  instability  results.  Note  that  for  the  case  of  a  deforming  bubble,  Kang  and  Leal  (1987) 
obtain  0.9<  Wccr  <1  -0  while  our  results  appear  to  indicate  1. 1<  Wscr  <\.2.  This  discrepancy 
in  the  critical  value  of  Weber  numbers  is  somewhat  puzzling.  However,  as  pointed  out  in 
Ryskin  and  Leal  (1984),  the  critical  Weber  numbers  for  bubbles  {fi  0,X  0)  and  drops 
are  expected  to  be  different.  This  is  due  to  the  non-negligible  density  and  viscosity  of  the 
fluid  in  the  drops.  In  our  computations,  we  expand  the  parameter  range  of  density  and  viscos- 
ity combinations.  To  gain  more  insight,  the  following  description  is  given. 

Considering  the  normal  stress  balance  at  the  interface  represented  by  Eq.  (2.12),  the 
pressure  on  each  side  of  the  interface  needs  to  be  obtained.  This  is  achieved  if  one  defines 
a  system  of  coordinates  at  the  interface,  such  that  r]  and  ^  are  coordinates  oriented  along  the 
tangential  and  normal  directions  to  the  interface.  Now,  the  pressure  in  each  phase  at  the  inter- 
face may  be  obtained  by  solving  the  momentum  equation  at  the  interface,  which  may  be  writ- 
ten in  the  form  (Ryskin  and  Leal  1984) 
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Figure  3.11.  Deformation  of  a  bubble  under  uniaxial  extensional  flow  at  Re=  1 0. 
(a)  Curve  shows  the  deformation  parameter  D  vs  Weber  number  We.  The 
critical  Weber  number  for  this  case  lies  between  1 . 1  and  1 .2.  (b)  L  vs  t  plot  for 
three  different  Weber  numbers  (We=0.7,  1,1,  1.2). 


62 


V{  p  +  qVV  )  ={QVX(D)-2Vxaj 


(3.9) 


where  (o  is  the  vorticity.  Integrating  this  expression  along  the  T^-direction  and  considering 
the  idealized  case  of  inviscid,  irrotational  flow,  one  obtains: 


where  the  subscript  i/o  implies  that  the  expression  applies  to  both  the  inner  and  outer  phases 
of  the  interface  and  Vfj  is  the  rj  (tangential)-component  of  the  velocity  at  the  interface.  The 
component  of  pressure  Pc  is  hydrostatic.  The  difference  between  the  inner  and  outer  hydro- 
static pressure  components  is  adjusted  as  described  previously  for  enforcing  mass  conserva- 
tion. Under  this  idealization,  the  dynamic  pressure  plays  the  dominant  role  in  the  normal 
stress  balance  at  the  interface.  This  is  an  important  term  that  distinguishes  the  behavior  of 
bubbles  from  the  inviscid  drops.  In  the  case  of  the  bubble,  the  dynamic  contribution  to  the 
pressure  from  the  inner  phase  is  negligible  and  the  stability  of  the  bubble  is  influenced  only 
by  the  outer  dynamic  contribution.  For  the  case  of  an  inviscid  drop  with  the  densities  ^=1, 
the  kinematic  condition  implies  fV^j,  =  (V^jo.  Therefore,  under  this  condition,  when  the 
pressures  at  the  interface  as  given  by  Eq.  (3.10)  are  inserted  into  Eq.  (2. 12),  the  dynamic  con- 
tributions to  the  normal  stress  balance  exactly  cancel.  In  inertia  dominated  flows  it  is  the 
dynamic  contribution  that  dominates  in  the  stress  balance.  Therefore,  in  this  situation  the 
behavior  of  the  inviscid  drop  will  be  vastly  different  from  that  of  the  bubble  where  the  inertia 
effect  is  negligible,  and  different  densities  of  the  drop  will  cause  different  critical  Weber 
numbers. 

For  the  case  of  a  drop,  since  the  dynamic  contributions  from  the  inner  surface  partially 
counteracts  the  dynamic  pressure  contribution  on  the  outer  phase,  the  critical  Weber  number 
may  be  expected  to  be  higher  than  that  for  a  bubble.  In  the  simulations  presented  here,  the 
ratios  for  density  and  viscosity  are/3=0.02  andX=0.02.  To  see  what  the  flow  field  character- 
istics are  in  this  case,  we  plot  in  Figure  3.12(a)  the  velocity  vectors  for  a  We=L2,  i.e.,  the 
critical  value.  The  pressure  contours  and  profiles  are  shown  in  Figures  3. 12(b)  and  (c),  re- 


Pi/o  = 


(3.10) 
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Figure  3.12.  flow  field  characteristics  for  the  case  of  drop  in  extensional  flow 
at  Re=10,  We=1.2,  |3=0.02,  X,=0.02.  (a)  Velocity  vectors  in  the  region 
containing  the  drop,  (b)  Pressure  contours  in  the  same  region,  (c)  Pressure 
profiles  at  8  different  sections  along  the  x-axis  starting  from  i=4  to  i=95  (drop 
extends  from  i=l  toi=98).  Arrow  indicates  increasing  i.  (d)  u-velocity  profiles 
at  the  same  stations  as  in  (c). 
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spectively,  at  8  equally-spaced  sections,  starting  from  i=4  to  95,  along  the  x-axis  in  the  re- 
gion occupied  by  the  drop  that  extends  from  i=l  to  98.  The  arrow  in  Figure  3. 12(c)  shows 
the  direction  of  increasing  /.  As  seen  in  Figures  3. 12(b)  and  (c),  there  is  ajump  in  the  pressure 
at  the  interface  location.  It  is  seen  that  although  the  pressure  at  each  x-station  inside  the  drop 
is  practically  constant,  there  is  a  variation  of  the  pressure  in  the  drop  in  the  tangential  direc- 
tion. In  fact,  the  pressure  rises  somewhat  as  one  goes  from  the  top  of  the  drop  toward  the 
side.  Although  this  rise  in  pressure  is  not  as  high  as  that  on  the  outer  phase  as  can  be  clearly 
shown  in  Figure  3.12(c),  this  may  lead  to  the  enhanced  stability  in  the  present  case  when 
compared  to  a  bubble  where  the  pressure  inside  is  a  constant  and  is  solely  hydrostatic.  Note 
that  in  the  interior  of  the  drop  the  pressure  along  the  interface  rises  from  the  value  of  about  1 
1.7  to  1.95  {dp=0.25),  while  outside  the  drop  it  decreases  from -0.8  to  about -2.3  (dp=-l.5). 
Therefore,  the  tangential  increment  of  the  pressure  inside  the  drop  is  non-negligible 
compared  to  that  outside.  , 

In  Figure  3. 12(d)  we  present  the  M-velocity  profiles  at  the  same  sections  as  the  p-pro- 
files  in  Figure  3. 12(c).  The  arrow  in  Figure  3. 12(d)  also  indicates  the  direction  of  increasing 
/.  As  can  be  seen  large  velocity  gradients  exist  in  the  interior  of  the  drop  due  to  the  strongly  I 
recirculating  flow  present  there.  However,  the  viscous  stresses  in  the  drop,  in  relation  to  that 
outside,  may  not  significantly  influence  the  drop  stability  via  Eq.  (3.12)  in  view  of  the  iner- 
tia-dominated flow  and  due  to  the  viscosity  inside  the  drop  being  0.02  times  the  outer  viscos- 
ity. It  is  not  straightforward  to  quantify  the  individual  contributions  of  the  viscous  stresses 
and  pressure  terms  in  the  balance  expression  given  by  Eq.  (2. 12)  due  to  the  non-linearity  of 
the  flow  field.  At  the  present  time,  based  on  the  flow  quantities  displayed  in  these  figures, 
it  is  proposed  that  the  discrepancy  in  the  Wccr  value  is  due  to  the  fact  that  even  for  the  fairly 
large  density  and  viscosity  jumps  considered  here,  the  dynamics  internal  to  the  drop  is  non- 
negligible.  Further  investigation  of  these  aspects  is  required  and  is  in  progress.  Suffice  it  I 
to  say  that,  as  pointed  out  by  Ryskin  and  Leal  (1984),  the  results  from  studies  of  the  behavior 
of  bubbles  do  not  carry  over  directly  to  inviscid  drops,  except  at  negligible  Reynolds  num- 
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bers;  for  the  imposed  flow  at  small  Reynolds  numbers,  the  pressure  inside  the  bubble  and 
inviscid  drop  are  constant,  since  there  is  no  dynamic  contribution  to  pressure. 

3.2.3.  Deformation  of  Drops  in  Constricted  Tubes 

In  the  above  sections  we  have  demonstrated  the  ability  of  the  algorithm  to  1)  handle 
complex  flows  in  the  arbitrary  domains  represented  on  a  Cartesian  grid,  2)  to  accurately  ob- 
tain the  dynamics  of  interfaces  under  the  balance  of  inertia,  viscous  and  capillary  forces.  In 
this  section  we  tie  these  aspects  together  in  studying  the  deformation  of  drops  in  viscosity 
and  inertia  dominated  flows  through  a  constricted  tube. 

The  interaction  of  drops  with  the  geometry  in  which  they  are  constrained  to  traverse  has 
received  some  attention  from  experimental  (Olbricht  and  Kung,  1992;  Olbricht  and  Leal 
1993)  and  numerical  viewpoints  (Tsai  and  Miksis,  1994;  Manga,  1996a,  b).  Since  the  major 
focus  has  been  on  the  behavior  of  multiphase  flow  through  porous  media  with  applications 
in  oil  recovery,  such  flows  have  been  restricted  to  the  creeping  flow  regime.  The  boundary 
integral  technique  is  applied  to  simulate  such  a  problem  (Tsai  and  Miksis,  1994).  In  recent 
work.  Manga  ( 1 996a,  b)  applied  the  same  technique  to  study  the  behavior  of  drops  in  geome- 
tries such  as  a  driven  cavity  and  bifurcating  channel,  again  restricted  to  Stokes  flow. 

To  our  knowledge,  numerical  work  on  the  behavior  of  drops  in  fixed  geometries  with 
significant  inertia  effects  has  been  restricted  to  the  study  of  drops  in  straight  capillaries  (Lee 
et  al.,  1995).  While  the  boundary  integral  technique  is  unable  to  handle  flows  with  convec- 
tion, purely  Eulerian  methods  have  to  contend  with  limitations  imposed  by  the  Cartesian  grid 
on  which  they  are  based.  The  present  algorithm  does  not  suffer  from  these  limitations  and 
thus  will  be  employed  below  to  investigate  the  flow  of  drops  in  a  constricted  tube  at  non-van- 
ishing Reynolds  numbers. 

Figure  3.13  shows  the  results  of  the  simulation  of  drop  motion  through  a  constricted 
tube  at  small  Reynolds  number.  For  this  case,  we  have  employed  the  same  conditions  as  one 
of  the  cases  in  Tsai  and  Miksis  ( 1 994) .  The  constriction  in  the  tube  is  a  sine-wave  with  ampli- 
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(a) 


Figure  3.13.  Flow  of  a  drop  through  a  constricted  tube  for  the  case  of  Re=0.0, 
Ca=0.1,  (3=1,  X=0.01,  a=0.9.  (a)  The  interface  shapes  shown  at  non-dimensional 
times  of  0.0,  0.25,  0.5,  0.75,  1.0,  1.25,  1.5,  1.75,  2.0.  (b)  ,(c)  and  (d)  Velocity 
vectors  for  the  times  0.25,  1.25  and  1.5  respectively. 
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tude  0.6R  and  wavelength  2.0  units,  where  R  is  the  radius  of  the  tube.  For  this  case  Ca=0.1, 
the  density  ratio  ^  =  Qi/Qo  =  i,  the  viscosity  ratio  X  =  jUi/jUo  =  0.01,  the  ratio  of  the  unde- 
formed  drop  radius  to  the  tube  radius,  a=0.9.  As  in  the  case  of  Tsai  and  Miksis  (1994),  the 
circular  drop  is  placed  initially  such  that  its  center  is  at  x=4.0.  The  initial  shape  of  the  drop 
is  shown  by  the  dotted  lines  in  Figure  3. 13(a).  The  evolution  of  the  drop  is  shown  at  equal 
intervals  (<5^=0.25)  of  time.  The  interface  shapes  at  these  instants  agree  remarkably  well  with 
the  results  of  Tsai  and  Miksis  (1994)  using  the  boundary  integral  technique.  As  can  be  seen, 
the  drop  deforms  essentially  in  conformity  with  the  shape  of  the  constriction  for  this  case. 

In  Figures  3.13(b)  to  (d)  we  show  the  velocity  vectors  for  this  case  at  instants  before, 
during  and  after  passage  of  the  drop  through  the  constriction.  The  effect  of  capillarity  can 
be  seen  in  these  figures  from  the  turning  of  the  velocity  vectors  in  the  region  where  the  flow 
field  is  not  strong,  such  as  at  the  regions  close  to  the  wall.  Also  to  be  noticed  is  that  at  the 
constriction,  the  front  of  the  bubble  is  elongated  whereas  the  rear  is  flattened  by  the  accelera- 
tion of  the  flow  at  that  location. 

As  an  additional  element  in  the  physics  regarding  to  how  the  inclusion  of  inertia  alters 
the  dynamics  of  the  drop,  we  simulate  the  flow  of  the  same  drop  in  the  constricted  tube,  but 
under  convective  flow  conditions.  If  the  convection  is  not  negligible,  the  boundary  integral 
technique  used  by  Tsai  and  Miksis  (1994)  cannot  be  applied  due  to  the  non-linear  field  equa- 
tions. The  Reynolds  number  based  on  the  tube  radius  is  50.  For  the  case  shown  in  Figure 
3.14,  Ca=0.1,p=l,X=0.01  and  a=0.9.  The  response  of  the  drop  in  this  situation  is  signifi- 
cantly different.  Due  to  the  presence  of  the  recirculation  zone  aft  of  the  constriction  the  pic- 
ture changes  drastically.  For  the  area  ratio  and  Reynolds  number  considered,  there  appears 
to  be  an  unsteady  motion  in  the  region  following  the  constriction.  A  region  of  recirculating 
flow  is  formed  and  travels  downstream  of  the  throat.  This  region  of  the  flow  has  a  profound 
effect  on  the  interface  shape.  The  interface  gets  stretched  into  an  elongated  shape  down- 
stream of  the  constriction  due  to  the  higher  flow  speeds  there  on  account  of  the  rapid  area 
change.  Since  the  material  boundary  is  convected  with  the  flow  the  elongated  drop  shape 


Figure  3. 14.  Entry  of  a  drop  into  a  constricted  tube  for  Re=50,  Ca=0. 1,  p=l, 
X=0.01,a=0.9.  (a)  drop  shapes  at  t=0, 0.25, 0.625, 0.875,  l.Oand  1.12.  (bHf) 
Velocity  vectors  at  t=0.25,  0.625,  0.875,  1.0  and  1.12,  respectively. 
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is  formed.  As  the  drop  develops  however,  the  shedding  recirculating  flow  region  has  an  ef- 
fect on  the  interface  and  draws  the  front  end  of  the  interface  into  the  shape  seen  in  Figure 
3.14(a).  In  this  case,  as  the  interface  deforms,  the  rear  of  the  drop  is  pushed  toward  the 
constriction  and  flattened.  Ultimately,  the  drop  collides  with  the  wall,  at  which  point  we  ter- 
minate the  calculations  due  to  uncertainty  in  the  subsequent  behavior. 

One  possibility  is  that  the  drop  will  eventually  pass  through  the  constriction,  as  the  film 
of  liquid  squeezed  between  the  drop  and  the  wall  of  the  tube  may  act  as  a  lubricating  layer, 
facilitating  the  motion  of  the  interface  relative  to  the  wall.  This  behavior  cannot  be  captured 
due  to  the  resolution  used  by  the  grid.  Alternatively,  the  collision  with  the  wall  may  result 
in  rupture  of  the  drop  surface  leading  to  breakup.  Since  the  impending  events  are  somewhat 
nebulous,  we  present  results  only  up  to  the  stage  where  the  collision  is  detected  by  the  algo- 
rithm. 

As  a  contrasting  case,  we  perform  the  flow  simulation  with  convective  effects  on  a 
smaller  drop,  for  which  a=0. 75.  All  other  conditions  are  identical  to  the  previous  case.  From 
Figure  3. 15,  it  is  seen  that  the  dynamics  is  very  similar  to  the  case  of  the  larger  drop,  except 
that  eventually  the  drop  passes  through  the  constriction.  Since  the  calculation  was  carried 
further  than  in  the  previous  case,  it  is  seen  that  the  translating  recirculating  flow  structure 
downstream  of  the  constriction  continues  to  deform  the  drop  resulting  in  the  flattened  struc- 
ture of  the  fore  end  of  the  drop.  For  the  given  Ca=0.1,  the  inertia  effects  appear  to  dominate 
over  the  capillary  forces  as  there  is  hardly  any  difference  in  the  behavior  of  the  smaller  drop 
in  comparison  to  the  larger  drop.  Due  to  the  larger  surface  curvatures  of  the  smaller  drop, 
in  the  case  of  sufficiently  small  Capillary  numbers,  it  is  expected  that  the  drop  size  will  have 
an  effect  on  its  interaction  with  the  flow  field.  Also,  due  to  the  unsteady  nature  of  the  flow 
field,  the  precise  shape  of  the  drop  downstream  of  the  throat  will  depend  on  the  relative  phase 
of  shedding  in  that  region  and  drop  entry  aft  of  the  throat.  Also,  from  preliminary  studies 
for  a  lower  Reynolds  number,  where  shedding  is  absent,  the  behavior  of  the  drop  is  noticed 
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Figure  3.15.  Entry  of  a  drop  into  a  constricted  tube  for  Re=50,  Ca=0.1,  P=l,  >t=0.01, 
a=0.75.  (a)  Shapes  of  the  drop  at  t=0.0,  0.25,  0.5,  0.625,  0.875,  1.25,  1.375,  1.5. 
(b)-(g)  Velocity  vectors  at  t=0.25,  0.5,  0.625,  0.875,  1.25,  1.5,  respectively. 


to  be  significantly  different  from  that  presented  here.  These  studies  are  currently  being  in- 
vestigated but  are  not  part  of  the  scope  of  the  present  thesis. 

3.3.  Summary 

In  this  chapter,  we  have  focused  on  the  numerical  aspects  of  the  algorithm  and  have 
tested  several  well-defined  and  challenging  flow  problems.  The  performance  of  the  method 
has  been  demonstrated.  The  applications  include: 

( 1 )  .  Flows  in  complex  2-dimensional  geometries  under  significant  convection,  which  obvi- 

ates the  need  for  grid-generation.  In  this  section,  we  have  compared  results  where  pos- 
sible to  those  given  by  a  boundary-fitted  code  and  shown  good  comparison. 

(2)  .  The  results  of  the  deformation  and  recovery  of  viscous  drops,  subjected  to  an  uniaxial 

extensional  flow,  show  good  agreements  with  those  of  Stone  and  Leal  (1989a,  b).  Thus, 
it  is  established  that  the  viscous-capillary  balance  is  well  captured  by  the  numerical 
technique. 

(3)  .  The  deformation  of  drops  under  uniaxial  extension  flow  under  convective  conditions 

has  been  studied.  In  contrast  with  the  work  by  Kang  and  Leal  (1987),  here  we  consider 
the  motion  of  the  fluid  inside  the  drop.  It  is  interesting  to  note  that  the  estimated  critical 
Weber  number  for  the  present  simulation  is  somewhat  higher  than  that  predicted  by 
Kang  and  Leal  (1987)  for  a  bubble  with  no  internal  dynamics.  As  discussed  in  Reskin 
and  Leal  (1984),  this  discrepancy  in  predicting  critical  Weber  numbers  is  expected.  In 
line  with  Mashayek  and  Ashgriz  (1996),  the  relaxation  oscillations  of  the  drop  under 
the  competition  between  inertia  and  capillary  forces  is  observed. 

(4)  .  Finally,  the  various  aspects  of  the  algorithm  tested  above  are  grouped  together  in  order 

to  study  the  behavior  of  drops  in  constricted  tubes.  For  low  Reynolds  number  flows, 
we  show  that  the  results  of  Tsai  and  Miksis  (1994)  are  exactly  reproduced.  But  we  go 
beyond  this  creeping  flow  regime,  and  find  some  interesting  behaviors  of  the  drop  in 
the  region  downstream  to  the  throat  at  Re=50. 


CHAPTER  4 
DYNAMICS  OF  COMPOUND  DROPS 

The  focus  of  this  chapter  is  to  analyze  the  deformation  and,  in  particular,  subsequent 
recovery  of  3-layer  or  compound  drops  in  order  to  shine  some  light  on  our  current  under- 
standing of  the  leukocyte  rheology. 

Recovery  experiments  have  been  extensively  used  in  recent  years  in  order  to  character- 
ize the  rheological  properties  of  leukocytes  (white  blood  cells).  In  these  experiments,  a  cell 
of  about  8  [im  in  diameter  is  aspirated  into  and  expelled  from  a  3  ~  7  [xm  diameter  glass  pipet. 
The  advantages  of  the  recovery  experiments  over  those  involving  the  flow  of  a  cell  into  a 
pipet  are  that  they  are  easier  to  perform  and  analyze,  and  that  the  potential  adhesion  of  the 
cell  to  the  pipet  is  not  an  issue.  Using  micropipet  techniques,  mechanical  properties  of  leuko- 
cytes have  been  evaluated  with  Newtonian  and  non-Newtonian  drop  models.  However,  re- 
ported experiments  indicate  that  these  existing  models  fail  to  explain  the  complex  rheologi- 
cal behavior  of  leukocytes  because  they  do  not  consider  the  internal  structure  of  the  cells. 
It  has  been  suggested  that  a  more  adequate  model  for  describing  the  cell's  behavior  is  the 
compound  liquid  drop  model  (Hochmuth  et  al.,  1993). 

4. 1  ■  Problem  Statement  and  Mathematical  Formulation 

We  consider  the  dynamics  of  a  compound  drop  immersed  in  a  uniaxial  extensional  flow 
and  its  subsequent  recovery  when  the  imposed  flow  is  removed.  The  physical  situation  cor- 
responds to  three  incompressible  Newtonian  fluid  layers  of  density  Qj  and  viscosity  which 
occupy  the  regions  Qj  (i=l,  2,  3).  The  regions  Qi  ,  Q2  ^nd  Q3  represent,  respectively,  the 
suspending  fluid,  shell  and  core  as  shown  in  Figure  4. 1(a).  The  surface  tension  coefficients 
a\2  and  023  at  the  two  interfaces  and  r23  are  assumed  to  be  constants.  The  dimensions 
R2  and  R3  are  the  undeformed  radii  of  the  outer  and  inner  drops,  respectively. 
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Figure  4.1.  Illustration  of  the  parameters  in  the  computation  of  the 
compound  drop  deformation  and  recovery,  (a)  The  definition  of  model 
parameters,  (b)  Computational  set-up  and  boundary  conditions. 

The  choices  of  characteristic  velocity,  length  and  pressure  are  as  follows: 

u  =  g,  ^  =  R,  p  =  e 

where  p,  =  Hi  is  the  viscosity  of  the  suspending  fluid,  6  =  a\2  is  the  surface  tension  of  the 
cytoplasm-suspending  fluid  interface.  Then,  the  dimensionless  governing  equations  for  axi- 
symmetric  creeping  flows  are: 
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(1)  Continuity  equations  in  region  Qj  (i=l,  2,  3): 


0 


(4.1) 


(2)  Momentum  equation  for  suspending  fluid  (region  Qi): 


0  =  -  Vp  +  V^V  +  f 


12 


(4.2) 


and  for  each  of  the  enclosed  fluids,  namely  the  shell  (^22)" 


0  =    -  Vp  +  XgV^V  +  f,2  +  f; 


23 


(4.3) 


and  the  core  (Q3): 


0  =    -  Vp  +  XgV^V  +  f; 


23 


(4.4) 


In  the  above  equations,  X2  =  1^2/^1  is  the  viscosity  ratio  of  shell  fluid  to  suspending  fluid, 
X,3  =  ^3/(11  is  the  ratio  of  core  viscosity  to  suspending  fluid  viscosity,  Vj  and  p,  represent  the 
velocity  and  pressure,  respectively,  in  region  Qj.  The  forces  f\2  and  f23  denote  the  Eulerian 
contributions  of  the  normal  stress  discontinuity  due  to  the  capillary  forces  evaluated  on  the 
interfaces  and  r23.  These  forces  are  transmitted  from  the  interfaces  to  the  fluid  by  use 
of  the  Immersed  Boundary  Technique  (Peskin,  1977;  Unverdi  and  Tryggvason,  1992).  The 
form  of  the  force  contribution  is  given  by: 


where  nj  is  the  unit  normal  vector,  defined  as  positive  when  pointing  outward  from  each  re- 
gion i,  as  shown  in  Figure  4. 1(a),  Xj  =  V  •  iij  represents  the  local  mean  curvature  of  the  in- 
terface; X  and  are  the  location  of  points  in  the  flow  field  and  on  the  interface,  6  denotes 
the  delta  function,  72  =  O12/0  and  73  =  023/0  are  the  dimensionless  surface  tensions  for  the 
shell  and  the  core,  respectively.  The  detailed  description  of  the  numerical  implementation 
of  these  force  terms  is  given  in  Chapter  2.4.2. 


12  ' 


23 


(4.5) 
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For  the  imposed  axisymmetric  uniaxial  extensional  flow,  the  far  field  condition  is 


Voo  =  E    R  =  ^ 


'2 
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0~ 
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0 

-  1 

y 

(4.6) 


where  E  is  the  strain-rate  tensor,  R  is  the  position  vector  and  G  represents  the  shear  rate  of 
the  imposed  flow.  Due  to  the  choice  of  the  velocity  scale,  the  nondimensional  boundary 
condition  at  infinity  becomes 


V    =  Ca 

Voo  2 
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(4.7) 


where  the  capillary  number  Ca  provides  a  measure  of  the  relative  dominance  of  viscous 
forces  to  surface  tensions  and  is  given  by  (fiiGR2/Oi2)- 

The  boundary  conditions  to  be  satisfied  on  each  interface  for  cases  involving  no  mass 
exchange  across  it  are  the  continuity  condition  and  the  stress  balances: 

(V)r„  =  (V),  =  (V)2 


(V)r,3  =  (V)2  =  (V)3 


(4.8) 
(4.9) 
(4.10) 
(4.11) 


"2  ■  T]  -  >^2"2  ■  ^2  =  Y2"2('^r,2  '  "2) 

"3  •  T2  -        •  T3  =  Y3"3(^r23  ■  "3) 

where  T  is  the  stress  tensor,  defined  as  follows. 

T=  -pI  +  [(VV)  +  (VV)T]  (4^2) 

It  is  noted  that  Eqs.  (4. 10-4. 11)  are  implicit  in  the  immerse  boundary  representation  of 
the  capillary  force  terms,  f  12  and  f23,  in  Eqs.  (4.2^.4).  The  advection  of  the  interface  is  per- 
formed with  the  kinematic  condition 


f  =  (Wr,  ■  Hi)  ~i 
where  X  represents  a  marker  point  on  the  interface. 


(4.13) 
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Considerations  regarding  accuracy  and  conservation  of  mass  in  each  region  of  the  com- 
pound drop  are  addressed  in  Chapter  2.5. 

4.1.1.  Computational  Setup 

The  choice  of  an  extensional  flow  to  pursue  our  study  is  motivated  by  the  need  to 
employ  a  simple  linear  flow  field  which  does  not  add  complexity  to  the  physics.  In  addition, 
the  behavior  of  single  phase  drops  (henceforth  simple  drops)  have  been  extensively  studied 
from  the  point  of  view  of  deformation  as  well  as  recovery  under  imposed  extensional  flows 
(Stone  et  al.,  1986;  Stone  and  Leal,  1989a,  b),  and  thus  provide  results  for  comparison.  Lin- 
ear uniaxial  extensional  flow  is  distinguished  from  simple  shear  by  the  absence  of  vorticity. 
No  net  translation  of  the  center  of  mass  of  the  drop  occurs  and  the  symmetries  dictate  that 
no  tank-treading  type  motion  is  possible. 

The  effects  of  this  straining  flow  on  a  compound  drop  are  investigated  by  solving  the 
Stokes  equations.  The  flow  conditions  and  computational  setup  are  indicated  in  Figure 
4.1(b).  The  symmetries  of  the  flow  are  exploited  to  mitigate  the  computational  effort  and 
only  the  first  quadrant  is  treated.  The  distribution  of  markers  along  the  interface  is  main- 
tained at  a  constant  arclength  spacing.  The  redistribution  along  the  interfacial  curve  is  per- 
formed at  every  time  step  to  prevent  undesirable  depletion  or  accumulation  of  markers.  The 
arclength  spacing  used  in  all  computations  presented  here  is  O.Bh,  which  was  found  by  exper- 
imentation to  provide  adequate  resolution  and  accuracy.  The  number  of  markers  on  the  inter- 
face varies  as  the  interface  stretches.  The  grid  distribution  is  nonuniform;  in  the  far  field, 
i.e.,  away  from  the  inner  region  shown,  the  grid  spacing  increases  linearly  toward  the  bound- 
aries of  the  domain,  while  in  the  inner  region  a  uniform  fine  grid  is  used.  This  yields  better 
resolution  in  the  desired  region  of  the  flow  field,  i.e.,  in  the  vicinity  of  the  interface.  The 
flow  field  is  imposed  as  shown,  thereby  exploiting  the  symmetry  of  the  problem.  The  incom- 
ing flow  at  the  top  and  outgoing  flow  at  the  right  are  held  constant  as  specified  by  Eq.  (4.7) 
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for  the  axisymmetric  case.  The  dimensions  L2  and  L3  are,  respectively,  defined  as  the  final 
deformed  extent  of  the  shell  and  core  in  the  x-direction  just  before  recovery  begins. 

4.1.2.  Selection  of  Model  Parameters 

Even  for  Newtonian  fluids,  the  compound  drop  model  is  represented  by  five  parameters 
(X,2>  ^3.  Y2>  Y3  and  R3)  as  opposed  to  two  (X,2  and  72)  for  a  simple  drop.  Apart  from  the  five 
model  constants,  an  additional  parameter  that  enters  in  this  study  is  the  capillary  number  of 
the  imposed  flow.  One  factor  that  can  be  fixed  is  the  value  of  R3.  This  is  based  on  exper- 
imental observation  that  the  nucleus  occupies  a  volume  of  2 1  %  of  the  neutrophil.  This  corre- 
sponds to  a  spherical  nucleus  of  nondimensional  radius,  R3=0.585.  Since  the  values  of  the 
other  parameters  and  the  nature  of  the  nucleus  and  shell  are  not  setded  (Dong  and  Skalak, 
1 992;  Hochmuth  et  al. ,  1 993 ;  Waugh  and  Tsai,  1 994),  we  assign  these  parameters  by  drawing 
on  experimental  observations.  It  is  generally  agreed  upon  that  the  nucleus  is  highly  viscous 
compared  to  the  shell.  Therefore,  we  will  examine  the  effect  of  introducing  a  viscous  core 
in  a  Newtonian  liquid  drop  and  identify  the  flow  and  deformation  conditions  that  are  mean- 
ingful to  explore. 

4.2.  Deformation  Analysis 

4.2. 1.  Steady-State  Configuration  (  Ca  <  Ca^-r ) 

We  first  examine  the  steady-state  shapes  of  the  simple  drop  against  the  results  of  Stone 
and  Leal  (1989a)  with  ^2=1  and  72=  1  •  In  the  case  of  a  simple  drop,  for  a  range  of  capillary 
numbers  below  the  critical  capillary  number,  Ca^r,  the  droplet  deformation  culminates  in  a 
steady-state  shape.  The  quantity  D=I(L-W)/(L+W)I  is  a  sensitive  indicator  of  the  final  drop 
shape  at  steady-state,  where  L  is  the  non-dimensional  major  radius  and  W  the  minor  radius 
of  the  ellipsoidal  drop  as  shown  in  Figure  4. 1(b).  This  quantity  therefore  provides  an  esti- 
mate of  the  aspect  ratio  of  the  steady-state  shape. 
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Figure  4.2.  The  steady-state  behavior  of  the  simple  and  compound  drops  when 
subject  to  a  uniaxial  extensional  flow,  (a)  D  vs.  Ca  for  the  simple  (^2=1,  Y2=l) 
compound  drops  (X,2=l,  X3=l,  Y2=1.  Y3=l>  R3=0.5):  •  :  simple  drop  (Stone  and  Leal, 

1989a);  :  simple  drop  (present  results);  o  :  compound  drop  (Stone  and  Leal, 

1990);  :  compound  drop  (present  results).  Shapes  corresponding  to  the  capillary 

numbers  I-»IV  are  shown,  (b)  Velocity  vectors  at  steady-state  for  Ca=0.1. 
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Our  results  for  a  simple  drop  with  nondimensional  viscosity  ^2=1  are  in  good  agree- 
ment with  those  of  Stone  and  Leal  (1989a),  including  the  prediction  of  the  critical  capillary 
number,  as  indicated  in  Figure  4.2(a).  This  latter  is  the  value  approached  asymptotically  by 
the  D  vs.  Ca  curves.  Beyond  the  critical  capillary  number,  no  steady-state  shape  is  possible 
and  the  drop  stretches  continuously.  The  critical  capillary  number  of  Cacr=0. 12  obtained  in 
the  present  case  corresponds  to  the  value  obtained  by  Stone  and  Leal  (1989a).  Also  shown 
in  Figure  4.2(a)  is  the  result  for  a  compound  drop  with  a  core  viscosity  >.3=1  and  radius 
R2=0.5,  and  a  shell  viscosity  X,2=l.  The  upper  solid  curve  corresponds  to  the  whole  com- 
pound drop,  while  the  lower  solid  curve  represents  the  core.  It  is  noted  that  although  the 
compound  drop  deforms  markedly  and  becomes  unstable  at  a  lower  capillary  number  than 
the  simple  drop,  the  deformation  of  the  core  is  very  limited.  Also  to  be  noted  is  that,  in  agree- 
ment with  Stone  and  Leal  (1990),  and  as  deducible  from  the  flow  field  shown  in  Figure 
4.2(b),  the  core  deforms  into  an  oblate  spheroid,  while  the  shell  becomes  prolate.  The  D  vs. 
Ca  curves  for  the  compound  drop  model  shows  behavior  in  excellent  agreement  with  those 
of  Stone  and  Leal  (1990).  The  steady-state  velocity  vectors  in  the  drop  are  shown  in  Figure 
4.2(b)  for  a  capillary  number  just  below  the  critical  value,  Ca=0. 1 .  The  recirculating  flows 
inside  the  compound  drop  can  be  seen  to  cause  the  oblate  and  prolate  deformation  of  the 
steady-state  inner  and  outer  interfaces. 

4.2.2.  Unsteady  Response  (  Ca  >  Catr ) 

In  this  section,  we  retain  the  condition  that  the  viscosity  of  the  shell  fluid  X,2=l,  while 
the  value  of  R2  changes  from  0.5  to  0.585  for  the  reason  described  in  Section  4. 1 .2,  and  ex- 
plore how  the  presence  of  a  core  alters  the  behavior  of  the  simple  drop. 

In  Figure  4.3(a),  we  show  the  curves  of  L  vs.  time  for  a  simple  drop,  a  compound  drop 
with  a  core  viscosity  ^3=1  and  a  compound  drop  with  a  highly  viscous  core  ^3=100.  The 
deformation  of  the  shell  under  the  imposed  extensional  flow  first  proceeds  rapidly  for  all 
three  cases.  Thereafter  the  deformation  slows  down  leading  to  a  linear  growth  region.  After 
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Figure  4.3.  Response  of  the  simple  and  compound  drops  to  the  uniaxial  extensional 
flow  for  Ca=0.14  >  Cacp  (a)  L  vs.  time,  (b)  Rate  of  deformation  dL2/dt. 
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a  certain  extent  of  deformation  there  is  a  second  region  of  accelerated  stretching.  When  the 
rate  of  deformation  of  the  shell  (dL/dt)  is  plotted  against  time,  as  shown  in  Figure  4.3(b), 
it  can  clearly  be  seen  that  there  is  a  minimum  in  the  deformation  rate.  The  slowdown  in  the 
intermediate  stage  of  the  drop  deformation  in  each  of  these  cases  is  due  to  the  formation  of 
large  curvatures  at  the  lateral  ends  of  the  drop.  Thus,  during  this  period  the  capillary  forces 
opposing  deformation  are  higher,  resulting  in  a  lower  rate  of  overall  deformation.  As  the 
shell  continues  to  deform  however,  the  equator  of  the  shell  gets  progressively  flatter  and  in- 
stability results  when  the  curvature  at  the  equator  begins  to  achieve  negative  values.  At  this 
point,  the  capillary  forces  assist  the  flow  field  in  deforming  the  drop,  at  least  in  the  regions 
close  to  the  equator.  Thus,  less  energy  needs  to  be  expended  by  the  flow  in  deforming  the 
interface  and  the  shell  is  stretched  continuously.  This  aspect  of  the  dynamics  has  been  dem- 
onstrated for  simple  drops  by  Stone  and  Leal  (1989a). 

The  only  distinction  between  simple  and  compound  drops  in  the  deformation  process 
is  the  time  scale  of  deformation.  The  inclusion  of  the  core  enhances  the  overall  stiffness  of 
the  compound  drop  as  can  be  seen  from  the  dL/dt  vs.  time  curves,  even  though  its  presence 
produces  a  continuous  stretching  mode  of  the  compound  drop  at  a  smaller  Cacr  than  for  the 
simple  drop. 

4.3.  Recovery  Analysis 

We  next  consider  the  recovery  characteristics  of  simple  and  compound  drops  for  two 
different  imposed  flow  strengths,  Ca=0. 14  and  Ca=0.35.  In  each  of  these  cases,  we  first  de- 
form the  shell  to  a  desired  stretched  state  by  imposing  an  extensional  flow  at  a  given  capillary 
number  and  then  switch  off  the  flow  to  allow  recovery. 

4.3.1.  Effect  of  the  Capillary  Number 

In  Figure  4.4,  we  show  the  deformed  and  recovering  shapes  of  the  compound  drop  for 
the  above  two  capillary  numbers.  The  case  of  a  less  viscous  shell  (X.2=l)  and  highly  viscous 
core  (X.3=  100)  is  studied.  The  deformed  shapes  under  the  imposed  flows  are  shown  in  Figure 
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4.4(al )  for  Ca=0. 14  and  Figure  4(b  1 )  for  Ca=0.35.  We  define  that  tj  is  the  time  for  deforma- 
tion to  the  desired  extent  and  tr  is  the  time  for  recovery  starting  from  the  initial  deformed 
shape.  It  is  noted  that  the  characteristic  times  for  deformation  are  substantially  different  for 
the  two  cases.  Furthermore,  although  the  extent  of  deformation  is  the  same  (L2=3),  the 
shapes  are  different.  For  the  larger  rate  of  stretching,  Ca=0.35,  the  shell  assumes  a  sharper 
shape  laterally  due  to  stronger  imposed  shear  rate  effect. 

The  time  scale  of  deformation  of  a  drop  is  proportional  to  the  drop  viscosity  for  highly 
viscous  drops,  as  can  be  expected  from  the  physical  consideration  that  the  larger  of  the  visco- 
sities between  the  drop  and  suspending  phase  should  influence  the  dynamics  (Rallison, 
1984).  This  scaling  arises  naturally  by  obtaining  the  integral  equation  form  for  the  interface 
velocity  as  in  the  boundary  integral  methods  (Rallison  and  Acrivos,  1978).  A  natural  param- 
eter that  emerges  as  a  convenient  scaling  factor  for  the  time  of  deformation  of  a  highly  vis- 
cous simple  drop  is  (1+X)  where  X  is  the  drop  viscosity  non-dimensionalized  with  respect 
to  the  suspending  phase  viscosity.  Thus,  for  the  compound  drop  configuration  shown  in  Fig- 
ure 4.4,  the  response  of  the  shell  to  deformation  is  much  more  rapid  than  that  of  the  core  (cor- 
responding to  X3/X,2=  100).  Therefore,  at  the  higher  stretch  rates,  the  deformation  of  the  drop 
is  governed  by  the  less  viscous  shell. 

It  is  important  to  remember  that  the  recovery  rate  of  the  drop  depends  on  the  history 
of  deformation  for  both  simple  and  compound  drop  models.  It  is  recalled  that  although  the 
flow  field  itself  is  quasi-stationary,  the  time-history  of  deformation  is  stored  in  the  drop  in 
terms  of  its  shape.  This  is  obvious  from  Figures  4.4(al)  and  (bl).  Once  the  flow  is  switched 
off  in  the  recovery  phase,  the  return  of  the  shell  to  its  spherical  unstressed  condition  is  con- 
trolled by  the  starting  shape  at  release.  In  the  case  of  the  shell  deformed  at  Ca=0.35,  the  pro- 
duced surface  curvatures  are  higher  than  those  of  the  case  at  Ca=0.14,  hence  the  capillary 
forces  causing  recovery  are  larger.  This  accounts  for  the  faster  recovery  of  the  shell  after 
high  rate  of  stretching. 
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(al)  td=35,  tr=0  (bl)  td=5.5,  t,=0 


Figure  4.4.  Effect  of  the  capillary  number  or  imposed  flow  strength  on  the 
deformation  (L2=3)  and  subsequent  recovery  of  the  compound  drop  (72=!'  Y3=l' 
^2=1,  X3=100).  (a)  Recovery  shapes  for  the  compound  drop  deformed  at  low  strain 
rate  (Ca=0. 14).  tj  is  the  time  for  deformation  to  the  extent  shown,  tr  is  the  recovery 
time  starting  from  the  shape  in  (al).  (b)  Recovery  shapes  for  the  drop  deformed  at 
higher  strain  rate  (Ca=0.35).  The  sequences  of  figures  in  (a)  and  (b)  correspond  to 
identical  initial  extension  lengths. 
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4.3.2.  Effect  of  the  Core  Viscosity 

The  behavior  of  a  compound  drop  is  now  examined  for  two  core  viscosity  values,  X3=l 
and  X3=100.  The  shell  and  suspending  fluid  retain  identical  viscosities,  X2=X,i=l.  We  study 
the  recovery  characteristics  of  compound  drops  after  they  have  been  deformed  to  the  same 
extent  under  identical  imposed  flow  (Ca=0. 14).  Shown  in  Figure  4.5  are  the  velocity  vectors 
for  the  two  cases. 

In  Figures  4.5(al)  and  (bl)  we  show  the  velocity  vectors  just  before  the  flow  is  switched 
off  and  the  recovery  process  begins.  It  is  clearly  seen  that  the  less  viscous  core  is  only  slightly 
more  deformed.  In  both  cases,  we  have  maintained  the  same  surface  tension  value  for  the 
core.  Since  the  radius  of  the  core  is  0.585  times  the  whole  drop,  the  capillary  forces  on  the 
core  are  stronger  than  on  the  outer  shell  surface.  Thus,  the  core  is  harder  to  deform.  Further- 
more, during  the  stretching  process,  the  core  lies  in  a  region  of  low  flow  strength.  These 
combined  factors  make  the  value  of  the  core  viscosity  plays  a  relatively  minor  role  in  the 
deformation  stage.  Consequently,  the  time  taken  for  the  compound  drop  to  deform  to  the 
extent  shown  in  Figures  4.5(al)  and  (bl),  namely  L2=3,  is  not  very  different  in  the  two  cases 
(t(i=31  for  ^3=1  and  td=35  for  X3=100,  respectively). 

It  is  observed  that  during  recovery,  the  time  of  recovery  for  the  more  viscous  core 
(tr=48. 1)  is  much  longer  than  that  of  the  less  viscous  case  (tr=24).  This  is  due  to  the  fact  that, 
although  in  the  deformation  process  the  core  lies  in  a  region  of  low  flow  strength,  during 
recovery  this  is  not  the  case,  as  can  be  seen  in  Figure  4.5(a2).  The  capillarity  driven  flow 
attempts  to  deform  the  core.  For  the  less  viscous  core,  the  recovery  flow  is  able  to  deform 
it  into  an  oblate  shape  as  shown  in  Figure  4.5(a3).  Notice  that  although  in  Figure  4.5(a2) 
the  flow  vectors  in  the  core  are  in  the  same  direction  as  the  flow  in  the  shell,  in  Figure  4.5(a3) 
the  outer  flow  sets  up  a  circulation  in  the  opposite  direction  in  the  core.  The  velocity  magni- 
tudes in  the  final  stage,  i.e..  Figures  4.5(a5)  and  (b5),  are  small,  and  the  shapes  are  very  close 
to  spherical. 
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Figure  4.5.  Dependence  of  the  recovery  velocity  field  on  the  viscosity  of  the  core. 
Velocity  vectors  for  the  ^3=1  and  ^3=100  deformation  are  shown  in  (al)  and  (bl).  The 
flow  fields  in  the  recovery  stages  are  shown  in  (a2)-(a5)  for  ^3=1  and  (b2)-(b5)  for 
X,3=100.  The  coressponding  velocity  at  different  recovery  times  along  the  x-direction  is 
plotted  against  x  in  Figures  (a6)  and  (b6)  respectively.  At  tr=0,  L2=3  and  L3=0.618. 
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The  effect  of  the  highly  viscous  core  is  clearly  seen  by  comparing  Figures  4.5(a2)-(a4) 
with  Figures  4.5(b2)-(b4).  Figures  4.5(a6)  and  (b6)  show  the  velocity  along  the  x-direction, 
i.e.,  the  normal  velocity  of  the  recovering  interfaces  along  the  symmetry  axis  y=0.  It  can 
be  observed  that  the  velocities  at  the  enclosed  drop  are  much  smaller  than  at  the  outer  inter- 
face. This  is  particularly  the  case  for  the  highly  viscous  inclusion,  which  acts  like  a  blockage 
to  the  flow  induced  by  the  capillary  forces  at  the  outer  shell  interface.  The  U(x,y=0)  vs.  x 
curves  also  give  a  picture  of  the  flow  strength  at  the  different  stages  of  the  recovery,  indicat- 
ing for  example  that  the  velocity  magnitudes  shown  in  Figures  4.5(a5)  and  (b5)  are  very 
small  and  the  recovery  is  practically  complete.  The  velocity  field  in  the  vicinity  of  the  highly 
viscous  core  is  very  weak,  causing  it  to  remain  undeformed  in  the  recovery  process.  This 
enhances  the  overall  viscosity  of  the  compound  drop  and  causes  the  large  difference  in  the 
recovery  rates  mentioned  above. 

4.3.3.  Comparison  of  Simple  and  Compound  Drops 

We  now  compare  the  recovery  results  of  simple  and  compound  drops  for  the  cases  pres- 
ented above  by  plotting  the  L  values  versus  time.  In  Figure  4.6(a)  we  show  the  L  vs.  tr  curves 
for  the  capillary  number  (deformation  rate)  of  0. 14  and  deformation  extent  L2=2.  The  recov- 
ery of  a  simple  drop,  a  compound  with  ^3=1  and  a  compound  drop  with  X3=100  are  shown. 
The  difference  between  the  recovery  processes  of  different  drops  is  not  noticeable  in  the  ini- 
tial stage  of  recovery.  After  this  stage  the  effect  of  the  core  becomes  more  noticeable  in  the 
recovery  curves.  Thus,  the  compound  does  recover  much  like  a  simple  drop  in  the  initial 
stage;  but  as  the  recovery  proceeds,  the  presence  of  the  core  enhances  the  "apparent  viscos- 
ity" of  the  compound  drop.  Also,  as  the  viscosity  of  the  core  increases,  the  departure  from 
the  Newtonian  simple  drop  behavior  becomes  more  pronounced. 

Figure  4.6(b)  shows  the  recovery  characteristics  after  the  drop  has  been  subjected  to  the 
larger  deformation  rate  Ca=0.35.  In  this  case,  the  behavior  of  a  compound  drop  with  an  ex- 
tremely viscous  core  (X3=1000)  is  also  examined  in  order  to  establish  a  trend.  The  results 
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Figure  4.6.  Recovery  curves  for  the  simple  and  compound  drops  after 
deformation  at  different  strengths  (capillary  numbers)  of  the  imposed 
extensional  flow.  Shown  here  are  L  vs.  time  curves  for  (a)  Ca=0.14  and  (b) 
Ca=0.35.  L2=2  for  both  cases  before  recovery  is  started. 
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for  the  two  higher  viscosities  are  very  close  and  such  highly  viscous  cores  can  be  taken  to 
approximate  a  solid  inclusion  within  the  compound  drop.  A  noteworthy  point  in  this  plot 
however  is  the  close  agreement  between  the  compound  drop  with  the  core  viscosity  of  ^3=1 
and  the  simple  drop.  As  can  be  seen  in  Figure  4.6(b),  the  core  of  the  compound  drop  with 
^3=1  is  highly  deformed  by  the  imposed  extensional  flow.  Thus,  in  the  recovery  process  the 
core  does  not  behave  like  a  passive  blockage  but  is  itself  recovering,  a  situation  no  different 
from  the  Newtonian  simple  drop.  The  only  noticeable  deviation  occurs  in  the  final  stages 
when,  as  seen  in  Figure  4.6(b),  the  core  has  been  driven  to  an  oblate  shape  and  has  to  recover 
its  spherical  shape.  In  that  stage,  the  dynamics  of  the  shell-core  combination  deviates  from 
a  Newtonian  simple  drop  behavior. 

Finally,  for  the  two  extents  of  drop  deformation  studied  (L2=2  and  L2=3)  at  Ca=0.14, 
the  recovery  curves  are  very  similar.  This  is  true  for  both  simple  and  compound  drops  with 
core  size  R3  <  0.6. 

4.4.  Summary 

From  the  above  study,  the  following  decisions  are  made  regarding  the  parameters  for 
the  3-layer  leukocyte  model  to  be  presented  in  the  next  Chapter: 

( 1 )  .  The  capillary  number  of  the  flow  employed  to  deform  the  compound  drop  does  not  have 

a  significant  impact  on  the  recovery  when  the  core  viscosity  is  sufficiently  high.  There- 
fore Ca=0.14  is  chosen  in  what  follows. 

(2)  .  The  core  viscosity,  when  high  enough,  does  not  change  the  dynamics  markedly,  as  seen 

from  Figure  4.6(b),  where  the  values  of  ^3=100  and  1000  show  little  difference  in  the 
recovery  curves. 

(3)  .  The  extent  of  the  drop  deformation  does  not  have  a  significant  effect  on  the  recovery 

if  the  compound  drop  is  sufficiently  deformed.  As  a  consequence,  an  extension  value 
of  L2=2. 1  will  be  selected  for  the  leukocyte  modeling  study  in  order  to  match  some  re- 
ported experimental  data. 


CHAPTER  5 
LEUKOCYTE  MODELING 

5.  L  Choice  of  Model  Parameters 

In  this  Chapter,  we  examine  the  effects  of  a  viscous  cytoplasm  enclosing  a  highly  vis- 
cous nucleus.  The  material  parameters  for  the  nucleus,  cytoplasm  and  suspending  fluid  as 
well  as  the  surface  tensions  at  the  inner  and  outer  interfaces  are  chosen  based  on  values  per- 
taining to  leukocytes. 

In  view  of  the  nebulous  picture  of  the  properties  of  the  nucleus  and  cytoplasm  of  leuko- 
cytes (only  apparent  values  being  known  based  on  global  properties  and  there  being  wide 
discrepancies  in  the  reported  values),  we  choose  to  select  only  representative  values  for  the 
material  properties  of  the  3-layer  model.  The  "apparent  viscosity"  of  the  leukocyte  as  mea- 
sured in  the  recovery  experiments  by  Tran-Son-Tay  et  al.  ( 1 99 1 , 1 994b)  is  around  800  Poise. 
The  suspending  phase  viscosity  is  usually  in  the  vicinity  of  0. 1  Poise.  Non-dimensionalizing 
viscosities  with  the  suspending  phase  value  gives  X,app=8000  for  the  overall  cell  viscosity. 

Numerically,  dealing  with  the  deformation  of  a  drop  with  such  a  large  viscosity  ratio 
is  extremely  taxing,  since  the  deformation  time  scale  is  very  large  (O(Xapp)).  In  order  to 
speed  up  the  calculations  we  set  the  suspending  phase  viscosity  to  1  Poise.  Thus  the  nondi- 
mensional  viscosity  Xapp  of  the  cell  (i.e.  whole  cell)  is  800.  This  is  justified  from  the  fact 
that  for  viscosity  values  high  enough,  as  will  be  shown  later,  there  is  no  significant  change 
in  the  dynamics  of  liquid  drops,  except  that  the  time  scale  of  deformation  is  scaled  in  propor- 
tion to  the  viscosity.  The  viscosity  of  the  cell  is  some  averaged  property,  based  on  the  viscos- 
ity of  its  contents.  Therefore,  with  a  21%  volume  for  the  nucleus  (corresponding  to  the  leu- 
kocyte nucleus)  and  assuming  that  the  nucleus  is  100  times  more  viscous  than  the  cytoplasm, 
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we  assign  viscosity  values  ^2=40  and  X3=4000.  In  order  to  reduce  the  parameter  space,  we 
assume  for  now  that  the  nondimensional  surface  tensions,  Y2=73=1.  The  flow  is  applied  at 
a  capillary  number  of  0. 14  and  shut  off  when  the  cell  reaches  a  deformed  length  of  L2=2. 1 . 

5.2.  Effect  of  a  Viscous  Cytoplasmic  Layer 

In  contrast  to  the  cell  investigated  in  Chapter  4.3.2,  where  the  cytoplasm  fluid  viscosity 
was  equal  to  that  of  the  suspending  fluid  (i.e.,  ^2=1)>  the  cytoplasmic  fluid  viscosity  under 
study  here  is  much  larger  than  the  suspending  phase  (i.e.,  X,2=40).  The  recovery  curves  of 
the  cytoplasm  and  slightly  deformed  nucleus  are  represented  in  Figure  5.1(a)  by  the  solid 
lines  on  the  L(t)  vs.  tr  plot.  The  recovery  process  is  considered  completed  when  the  cell 
length  L<  1.05. 

The  recovery  curve  for  a  simple  Newtonian  drop  with  a  viscosity  value  equals  to  the 
cell  volume  averaged  viscosity,  namely  832,  is  also  shown  for  comparison  in  Figure  5. 1(a). 
The  recovery  process  of  a  simple  drop  is  much  slower  than  that  of  a  cell  with  a  comparable 
averaged  viscosity  value.  Thus,  the  use  of  the  volume  averaging  procedure  is  highly  inade- 
quate for  obtaining  the  effective  or  "apparent"  viscosity  of  the  cell.  In  fact.  Figure  5.1(a) 
shows  that  the  recovery  behavior  of  a  simple  drop  with  a  lower  viscosity  value,  e.g.,  240  or 
100,  does  approach  better  that  of  the  cell  under  study. 

Although  theoretical  estimates  exist  for  the  effective  viscosity  of  a  dilute  suspension 
of  compound  drops  (Stone  and  Leal,  1990),  the  only  measure  of  the  effective  viscosity  of 
a  compound  drop  is  qualitative.  The  only  guidance  offered  (Johnson  and  Sadhal,  1985; 
Stone  and  Leal,  1990)  is  that  the  effective  viscosity  of  the  drop  is  intermediate  between  that 
of  a  solid  sphere  (which  limit  is  approached  as  1),  and  that  of  a  simple  drop  (composed 
only  of  the  cytoplasm,  this  limit  of  course  being  realized  as  R3-»0).  Furthermore,  the  effec- 
tive or  "apparent  viscosity"  of  a  compound  drop,  as  measured  by  the  rate  at  which  the  drop 
recovers,  depends  on  the  deformation  history,  in  addition  to  its  material  properties. 
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Figure  5.1.  The  recovery  behavior  of  the  3-layer  leukocyte  model  compared  to 
several  averaged  simple  drop  models,  (a)  L  vs.  unsealed  recovery  time  tf.  (b)  L 
vs.  scaled  recovery  time  t^.  (c)  Scaled  rate  of  deformation  dL/dts  vs.  scaled 
recovery  time  ts.  (d)  L  vs.  scaled  recovery  time  tg.  All  viscosities  and  scaling 
factors  chosen  for  3-layer  model  are  indicated  in  the  legends. 
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A  different  perspective  is  obtained  by  plotting  the  L(t)  curves  against  a  rescaled  dimen- 
sionless  time  t^  defined  as, 

1  +  >^app  (5.1) 

where  the  X^pp  is  either  the  apparent  viscosity  of  the  cell  or  the  viscosity,  X,  of  a  simple  drop.  , 
A  similar  representation  of  the  recovery  results  by  rescaling  the  recovery  time  is  also 
employed  by  Tran-Son-Tay  et  al.  ( 199 1 , 1994a).  In  that  work,  experimentally  obtained  re- 
covery data  of  leukocytes  are  compared  to  the  behavior  of  Newtonian  simple  drops  by  fitting 
the  numerically  generated  curve  to  the  experimental  data.  i 

Here,  for  the  3-layer  model,  the  scaling  factor  is  chosen  to  be  Xapp=285  based  on  match- 
ing asymptotic  behavior  of  the  cell  with  that  of  a  simple  drop.  This  yields  the  picture  in  Fig- 
ure 5. 1  (b).  The  simple  drop  curves  collapse  together  if  the  dimensionless  recovery  time  de- 
fined in  Eq.  (5. 1)  is  used.  The  result  is  expected  because  of  the  similarity  in  the  recovery  - 
paths  for  highly  viscous  drops.  The  only  driving  forces  in  the  recovery  process  are  the  capil- 
lary forces,  which  for  the  same  starting  shapes,  should  only  yield  differences  in  recovery 
times.  '• 

Further  insight  can  be  obtained  by  plotting  the  scaled  rate  of  recovery  dL/dts  against 
ts  as  in  Figure  5.1(c).  It  is  clearly  observed  that,  while  for  ts  <  around  0.2  the  recovery  of 
the  cell  is  extremely  rapid,  a  much  lower  recovery  rate  takes  over  at  later  times.  There  is  a 
distinct  change  in  the  recovery  behavior  of  the  cells.  Thus,  the  rapid  recoil  ascribed  to  leuko- 
cytes in  initial  recovery  stages  can  be  provided  by  such  a  3-layer  model.  This  rapid  phase 
is  of  course  due  to  the  fact  that  the  nearly  undeformed  stiff  nucleus  does  not  participate  in 
the  recovery  dynamics,  as  seen  in  Chapter  4.3.2.  Therefore,  the  initial  recovery  stage  is  dom- 
inated by  the  cytoplasm  region.  The  viscosity  of  this  region  is  low  and  hence  the  time  scale 
of  response  is  fast,  based  on  the  scaling  in  Eq.  (5.1). 

When  the  deformation  of  the  cytoplasmic  layer  decreases,  the  interaction  of  the  cyto- 
plasmic surface  with  the  nuclear  surface  becomes  more  important,  leading  to  the  dynamics 
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of  the  cytoplasm-nucleus  combination  which  makes  this  period  one  of  slow  recovery. 
Therefore,  even  within  the  Newtonian  fluid  consideration  there  is  a  possible  wide  disparity 
in  observed  cell  apparent  viscosity  depending  on  the  recovery  characteristics.  We  will  tackle 
this  issue  in  more  detail  later. 

If  the  scaling  factor  for  the  cell  is  chosen  to  be  ^app=40  (the  viscosity  of  the  cytoplasm) 
instead  of  285,  the  earlier  phase  of  the  recovery  will  be  matched  as  seen  in  Figure  5.1(d). 
As  seen  in  this  figure,  the  recovery  behavior  in  the  initial  stage  follows  the  slope  of  the  simple 
drops,  before  deviating  from  those  curves  at  about  ts=0.2.  Thus,  at  least  in  the  initial  stages, 
the  apparent  viscosity  of  the  cell  is  close  to  that  of  the  cytoplasm. 

5.3.  Effects  of  Nucleus  Deformation  and  Time  Scales 

As  mentioned  in  the  previous  section,  it  is  possible  to  capture  an  initial  rapid  recoil  of 
a  leukocyte  with  a  3-layer  Newtonian  model,  supporting  the  suggestion  of  Hochmuth  et  al. 
( 1 993).  This  is  seen  to  be  possible  due  to  the  preferential  flow  of  the  less  viscous  cytoplasmic 
liquid  in  the  recovery  process.  However,  one  important  fact  has  to  be  considered  at  this  junc- 
ture. Under  certain  conditions,  such  as  long  holding  times  or  low  deformation  rates,  the  re- 
covery of  the  leukocytes  ensues  like  a  Newtonian  simple  drop.  Based  on  our  findings  above, 
it  appears  that  the  existence  of  a  nucleus  leads  to  the  distinct  fast  and  slow  regimes  in  the 
recovery  process.  In  order  to  examine  the  conditions  under  which  a  cell  can  yield  a  Newto- 
nian simple  drop-type  recovery  behavior  and  to  clarify  the  essential  elements  of  the  required 
3-layer  model,  we  will  study  the  recovery  dynamics  of  two  sets  of  ratios  (72/^2)  (Y3/X3). 

5.3.1.  Incompatible  Nucleus  and  Cytoplasm  Time  Scales 

We  want  to  investigate  the  impact  of  the  deformation  of  the  nucleus  on  the  cell  recovery. 
One  difficulty  encountered  in  this  regard  is  that,  as  mentioned  before,  the  extensional  flow 
is  weak  in  the  vicinity  of  the  nucleus.  Also,  it  is  difficult  to  obtain  sufficient  deformation 
of  the  nucleus  due  to  its  inherent  stiffness,  on  account  of  the  higher  viscosity  and  smaller 
radius  amounting  to  larger  capillary  forces.  We  attempted  to  increase  the  strain  rate  on  the 
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cell  by  raising  the  Capillary  number  to  0.35.  This  failed  to  achieve  the  desired  deformation 
of  the  nucleus,  due  to  the  enhanced  rapidity  of  flow  of  the  cytoplasmic  fluid.  The  outer  inter- 
face deforms  rapidly  and  contacts  the  nucleus  at  which  time  the  calculations  are  terminated. 
The  nucleus  itself  is  slow  to  deform  and  remains  mildly  extended.  Therefore,  our  strategy 
to  test  the  effect  of  the  nucleus  deformation  is  to  start  the  recovery  process  from  given  initial 
shapes  of  the  cytoplasm  and  the  nucleus,  and  exploit  the  quasi-stationary  nature  of  the  flow 
field.  In  the  present  computations,  the  cell  has  values  72=  1 .  Y3=  ^  >  ^2=40,  and  X,3=4000.  This 
corresponds  to  a  highly  viscous  nucleus.  The  shape  of  the  drop  is  taken  to  be  that  obtained 
after  deformation  to  L2=2.1  under  Ca=0.14,  and  the  shape  of  the  nucleus  is  imposed. 

The  velocity  field  in  the  recovery  period  is  shown  in  Figure  5.2(a)  for  the  case  where 
the  nucleus  is  deformed  to  a  very  small  extent  (L3=0.618)  which  resulted  from  the  applica- 
tion of  the  extensional  flow  as  in  the  previous  section.  Clearly,  in  the  recovery  stage,  the 
response  of  the  cell  is  limited  to  the  outer  interface  only,  while  the  nucleus  itself  serves  mere- 
ly as  an  internal  solid-like  obstruction  to  the  flow  field.  It  is  this  characteristic  of  the  flow 
field  which  leads  to  the  rapid  recoil  response  of  the  cell.  The  flow  fields  for  the  larger  extent 
case,  L3=l,  are  shown  in  Figure  5.2(b).  The  deformed  nucleus  is  assumed  to  be  a  cylinder 
with  a  hemisphere  cap  at  both  ends.  One  quadrant  of  the  cell  is  illustrated  in  Figure  5.2(bl). 
The  dimensions  of  the  cylinder  and  hemisphere  are  assigned,  for  a  given  extent  of  deforma- 
tion, so  that  the  volume  of  the  nucleus  is  maintained  at  21%  of  the  cell.  The  initial  shape 
of  the  cell,  i.e.,  just  before  recovery,  is  the  one  given  in  Figure  5.2(al).  It  can  be  seen  that 
due  to  the  higher  viscosity  of  the  nucleus,  the  flow  resulting  from  the  recovery  of  the  outer 
interface  is  decelerated  as  it  approaches  the  nucleus,  leading  to  a  rapid  initial  recoil  phase. 

The  apparent  viscosity  values,  Xapp=285  for  L3=0.618  and  X,app=220  for  L3=l,  are  ob- 
tained by  using  the  scaling  in  Eq.  (5.1)  when  the  asymptotic  recovery  trend  of  the  cell 
matches  with  that  of  the  simple  drop  (Figure  5.3).  Figure  5.3  shows  the  recovery  curves  of 
the  simple  drop  and  the  cell  with  different  initial  nucleus  deformations.  The  off-Newtonian 
trends  of  the  cells  are  clearly  seen  when  the  nucleus  is  only  slightly  deformed.  In  the  later 
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Figure  5.2.  Velocity  vectors  at  different  scaled  times  tr  for  the  recovery  of  the 
3-layer  cell  (72=  1.  Y3=l>  ^2=40,  }i3=4000)  with  different  initial  nucleus 
deformations  and  the  same  overall  cell  deformation.  Scaled  times  at  which  the 
vector  are  shown  are  indicated  alongside  each  case,  (a)  For  L3=0.618,  L2=2. 1 . 
The  scaling  factor  Happ  for  time  is  285 .  (b)  For  L3=  1  and  L2=2. 1 .  The  scaling 
factor  Xann  for  time  is  220. 
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part  of  the  recovery  process,  cells  with  an  initially  deformed  nucleus  behave  as  simple  drops. 
This  can  be  seen  in  Figure  5.3  where  the  recovery  curve  for  L3=l  at  a  scaled  time  ts*=^0.7 
begins  to  approach  the  simple  drop  behavior. 

Figures  5.2  and  5.3  demonstrate  the  significant  role  the  extent  of  deformation  of  the 
nucleus  has  on  the  rheological  behavior  of  the  cell. 


0.5  1  1.5  2 

ts=t/(l+Xapp) 


2.5 


Figure  5.3.  Recovery  paths  for  cells  with  different  initial  nucleus  deformation 
and  a  highly  viscous  nucleus  (72=  1,  Y3=l'  ^2=40,  X3=4000). 


5.3.2.  Compatible  Nucleus  and  Cytoplasm  Time  Scales 

In  this  study,  we  seek  parameters  such  that  the  time  scales  of  the  cytoplasm,  Xcytopiasm 
=  Y2/^2.  and  the  nucleus,  Tnucleus  =  Y3/^3'  of  the  cell  are  compatible.  The  selected  material 
parameters  are  72=1,  Y3=10.  ^2=40,  and  X3=400.  Although  the  time  scale  compatibility 
holds  when  the  nucleus  and  cytoplasm  layers  are  considered  as  two  separate  drops  of  the 
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same  radius,  it  is  not  clear  what  the  relationship  will  be  when  the  two  drops  are  allowed  to 
interact  in  the  cell. 

Figure  5.4(a)  shows  the  recovery  flow  field  for  the  case  corresponding  to  L3=0.618. 
As  can  be  seen,  the  nucleus  remains  practically  undeformed  and  again  acts  as  an  obstruction 
as  the  cell  recovers.  The  flow  is  decelerated  in  the  vicinity  of  the  nucleus  due  to  the  increased 
resistance  from  the  nucleus.  In  contrast,  when  the  nucleus  is  stretched  to  L3=l,  although  the 
cell  is  deformed  to  the  same  extent  and  the  resulting  capillary-induced  flow  from  the  outer 
interface  is  of  the  same  magnitude,  the  recovery  flow  picture  is  changed  dramatically,  as  seen 
in  Figure  5.4(b).  Now,  as  seen  from  the  velocity  vectors,  the  flow  field  within  the  nucleus 
closely  follows  that  of  the  flow  induced  by  the  outer  interface.  It  is  observed  that  because 
the  strengths  of  the  recovery  flow  field  induced  by  the  outer  and  inner  interfaces  are  compa- 
rable, the  cell  recovers  almost  like  a  homogeneous,  Newtonian  simple  drop  when  L3=l. 
Thus,  the  response  of  the  cell  depends  not  only  on  the  material  properties  of  the  nucleus,  but 
also  its  extent  of  deformation. 

The  velocities  along  the  symmetry  line  (y=0)  during  the  recovery  process,  for  both  of 
the  L2  values  above  are  shown  in  Figure  5.5(a)  and  (b).  It  is  observed  that  because  the 
strengths  of  the  recovery  flow  field  induced  by  the  cytoplasm  and  nucleus  interfaces  are 
comparable,  the  cell  recovers  almost  like  a  homogeneous  liquid  drop  when  L2=1.0,  i.e.,  for 
the  large  nucleus  deformation.  On  the  other  hand,  for  the  slightly  deformed  nucleus 
(L3=0.618),  the  flow  is  decelerated  in  the  vicinity  of  the  nucleus  due  to  the  increased  resis- 
tance from  the  nucleus.  Thus,  although  the  cell  is  deformed  to  the  same  extent  and  the  result- 
ing capillary-induced  flow  from  the  outer  interface  is  of  the  same  magnitude,  the  response 
of  the  whole  cell  will  depend  on  the  nucleus  deformation.  The  response  of  the  nucleus  to 
the  incoming  recovery  flow  field  due  to  the  outer  interface  depends  not  only  on  the  material 
properties  of  the  nucleus,  but  also  on  its  extent  of  deformation. 

The  curves  in  Figure  5.6(a)  are  plotted  by  matching  the  asymptofic  behavior  of  the  cells 
in  each  case  to  that  of  the  simple  drop  and  obtaining  the  corresponding  apparent  viscosity 
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Figure  5.4.  Velocity  vectors  during  recovery  for  the  nucleus  with  a  time  scale 
compatible  with  that  of  the  cytoplasm  (72=!.  Y3=10'  ^2=40,  ^3=400)  for  different 
initial  shapes.  Scaled  times  t^  during  recovery  are  indicated  in  each  case,  (a)  For 
L3=0.618,  the  scaling  factor  Xapp  for  time  is  165.  (b)  For  L3=l,  the  scaling  factor 
Xann  for  time  is  70. 
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(a)  L3=0.618,  Y2=l,  Y3=10,  >^2=40,  ^3=400 
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(b)  L3=l,  Y2=l,  Y3=10,  >^2=40,  ^3=400 
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Figure  5.5.  The  velocity  profile  along  x-direction  U(x,y=0)  at  the  symmetry  axis  y=0. 
(a)  The  U  velocity  for  the  instants  shown  in  Figure  5.4(a).  The  velocities  of  the  nucleus 
and  that  induced  by  the  outer  layer  recovery  are  seen  to  be  incompatible  at  the  location 
of  the  nucleus.  The  nucleus  acts  like  a  simple  blockage  and  is  not  deformed 
significantly,  (b)  U(x,y=0)  for  the  case  in  Figure  5.4(b).  The  recovery  flow  induced 
by  the  outer  interface  is  well  matched  with  that  of  the  nucleus,  thus  ensuring  a 
compatible  recovery  path. 
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values.  These  are  Xapp=70, 117  and  165  for  L3=l,  0.75  and  0.618,  respectively.  In  contrast 
with  the  previous  case.  Figure  5.3,  where  the  recovery  time  scale  of  the  nucleus  is  incompat- 
ible with  that  of  the  cytoplasmic  layer,  the  trend  in  the  recovery  curves  is  not  apparent  here. 
In  particular,  the  L3=l  and  L3=0.75  curves  cross  over  the  L3=0.618  curve  in  the  later  stages 
of  recovery.  Notice  that  this  can  be  correlated  directly  with  the  time  when  the  nucleus  is  re- 
turning to  a  nearly  spherical  shape,  i.e.,  the  nucleus  recovery  has  begun  to  slow  down,  which 
decouples  the  outer  interface  recovery  from  that  of  the  nucleus.  The  curves  seem  to  indicate 
that  the  recovery  paths  are  distinctly  different  from  the  Newtonian  curve,  with  an  initial  rapid 
transient  and  final  slow  recovery  when  the  nucleus  returns  to  its  undeformed  shape. 

However,  an  entirely  different  viewpoint  emerges  when  the  initial  recovery  phase  is 
matched  to  that  of  the  simple  drop  in  each  case.  These  curves  are  shown  in  Figure  5.6(b). 
The  apparent  viscosities  obtained  from  the  scaling  are  this  time  )iapp=35, 42  and  47  for  L3=  1 , 
0.75  and  0.618,  respectively.  The  interesting  feature  to  note  is  that  the  recovery  curve  for  the 
case  with  the  most  deformed  nucleus  matches  quite  well  with  that  of  the  Newtonian  simple 
drop  for  a  considerable  extent  on  the  recovery  curve.  There  is  a  slow  final  transient  period 
when  the  cell  curve  deviates  slightly  from  the  Newtonian  behavior,  this  being  at  the  stage 
when  the  nucleus  has  recovered  almost  completely,  as  seen  from  the  lower  set  of  curves.  It 
is  also  noted  that  there  is  now  a  monotonic  trend  seen  in  the  curves,  i.e.,  when  the  nucleus 
is  less  deformed  the  cell  increasingly  deviates  from  a  Newtonian  behavior.  In  all  these  curves 
there  is  a  long  transient  part  that  corresponds  to  the  final  recovery  of  the  cell  when  the  nucleus 
has  reached  steady-state.  The  recovery  of  the  nucleus  takes  place  faster  due  to  the  larger 
capillary  restoring  forces  acting  on  it  on  account  of  its  larger  surface  curvature. 

From  the  apparent  viscosities  listed  above  for  the  curves  in  Figure  5.6(b),  it  is  somewhat 
puzzling  that  the  apparent  viscosity  of  the  cell  for  L3=l  is  lower  than  that  of  the  cytoplasmic 
component.  This  can  be  explained  by  the  fact  that  although  the  time  scales  for  recovery  of 
the  nucleus  and  the  cytoplasm  are  compatible,  the  surface  curvatures  of  the  nucleus  are  high- 
er. Thus,  the  recovery  rate  of  the  nucleus  is  faster  than  that  of  a  simple  drop  with  the  same 
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Figure  5.6.  Recovery  curves  for  the  case  (Y2=l,  73=10-  ^2=40,  X.3=400)  when 
nucleus  deformation  time  scale  is  compatible  with  the  cytoplasmic  layer  time  scale, 
(a)  L  vs.  scaled  time  ts,  scaling  factors  obtained  to  match  the  long-time  recovery 
asymptote  of  the  simple  drop  curve  are  shown  in  legend,  (b)  L  vs.  ts  scaled  to  match 
the  initial  recovery  stages,  obtained  with  scaling  factors  shown  in  the  legend. 
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cytoplasmic  viscosity  producing  this  remarkable  result.  Therefore,  the  3-layer  model  for 
sufficient  nucleus  deformation  can  produce  an  apparent  viscosity  that  is  lower  than  the  cyto- 
plasmic  viscosity. 

S  4  Discussion  And  Tmplir-^tinns  to  Leukocyte  Rheology 

It  is  clear  from  the  results  presented  above  that  the  rheological  behavior  of  leukocytes 
modeled  as  a  3-layer  structure  can  be  intimately  linked  to  the  internal  structure  and  material 
properties  of  the  cell.  The  proposed  model  can  explain  the  experimental  results  as  long  as 
the  cell  and  its  nucleus  have  compatible  recovery  time  scales.  This  condition  is  required 
since  the  cell's  response  is  observed  to  be  Newtonian  under  specific  situations.  Present  esti- 
mates of  the  properties  of  the  cell,  as  proposed  by  Skalak  et  al.  (1990)  place  the  values  of 
cytoplasmic  viscosity  at  approximately  200  Poise,  nuclear  viscosity  as  2000  Poise  the  elastic 
modulus  of  the  cortical  membrane  at  about  300  dyn/cm^  and  the  elastic  modulus  of  the  mem- 
brane complex  surrounding  the  nucleus  at  about  4000  dyn/cm2  (pong  et  al.,  1991).  Thus, 
one  obtains  the  ratio  Y3/Y2^10  and  }.3/^2^10,  which  indicates  compatible  recovery  time 
scales,  in  line  with  the  above  statement.  The  difficulty  in  obtaining  a  confirmation  of  the 
nucleus  deformation  characteristics  lies  in  the  inaccessibility  of  information  on  the  shape  as- 
sumed by  the  nucleus  in  the  micropipet  experiments. 

The  initial  rapid  entry  of  the  cell  into  micropipets  can  be  explained  in  the  3-layer  model 
by  the  fact  that  the  initial  rapid  response  represents  the  flow  of  the  cytoplasm  into  the  pipet 
while  the  less  deformable  nucleus  lags  behind.  It  has  been  observed  by  Zhelev  and  Hoch- 
muth  (1994)  that  the  maximum  apparent  viscosity  of  the  cell  during  entry  occurs  when  half 
the  nucleus  is  within  the  micropipet.  This  observation  supports  our  finding  that  the  deforma- 
tion of  the  nucleus  determine  the  cell  apparent  viscosity. 

Aspiration  experiments  also  indicate  that,  as  the  rate  of  suction  increases,  the  apparent 
viscosity  decreases,  which  has  been  speculated  to  be  due  to  a  shear-thinning  effect  (Need- 
ham  and  Hochmuth,  1990;  Tsai  et  al,  1993;  Waugh  and  Tsai,  1994).  From  the  3-layer  cell 
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point  of  view,  at  high  rate  of  deformations,  the  response  of  the  cytoplasm  is  rapid,  and  it  im- 
mediately responds  to  the  applied  forces  while  the  nucleus  deforms  little,  thus  leading  to  a 
lower  apparent  viscosity  (Figure  4.4).  In  our  computations,  this  is  observed  in  connection 
with  our  efforts  to  raise  the  capillary  number  of  the  imposed  extensional  flow  on  the  cell  from 
0. 14  to  0.35,  as  mentioned  in  Chapter  4.3. 1 .  These  flows  fail  to  deform  the  nucleus;  instead 
the  cytoplasm  flows  rapidly  until  the  gap  between  the  cell  boundary  and  the  nucleus  becomes 
very  small.         .  -  ,  . 

It  is  observed  that,  contrary  to  the  supposed  non-Newtonian  nature  of  the  cell,  the  vis- 
cosity of  the  flow  inside  the  micropipet  remains  the  same  independent  of  the  aspiration  pres- 
sure (Needham  and  Hochmuth,  1990).  In  terms  of  the  3-layer  model,  this  possibility  results 
from  the  fact  that  once  the  nucleus  and  overall  cell  are  deformed  to  their  equilibrium  configu- 
rations in  the  micropipet,  the  apparent  viscosity  of  this  3-layer  structure  remains  a  constant, 
since  the  apparent  viscosity  is  simply  a  function  of  the  cell  and  nucleus  shapes  and  their  mate- 
rial properties.  The  cell  apparent  viscosity  remains  constant  once  the  steady-state  shape  is 
assumed  in  the  micropipet  and  the  effect  of  the  aspiration  pressure  then  is  simply  to  influence 
the  flow  rate  of  the  compound  cell  body. 

When  the  cell  is  deformed  rapidly,  by  the  application  of  a  small  force,  the  cytoplasmic 
layer  of  the  cell  responds  immediately,  while  the  nucleus  remains  almost  undeformed.  Upon 
release  of  the  force,  the  cell  recovers  with  a  rapid  recoil  response  at  a  rate  corresponding  to 
the  apparent  viscosity  of  the  cytoplasmic  liquid.  This  fast  phase  is  then  followed  by  a  slower 
phase  when  the  cell  membrane  approaches  the  nucleus.  As  argued  by  Dong  et  al.  (1991), 
this  small  deformation,  rapid  recoil  behavior  should  yield  a  measure  of  the  cytoplasmic  vis- 
cosity. Based  on  their  model,  this  value  is  estimated  to  be  of  the  order  of 50-600  Poise.  How- 
ever, this  value  could  have  been  much  lower  and  more  physiological  if  it  was  determined 
with  the  3-layer  model.  A  correct  value  cannot  be  obtained  until  additional  information 
about  the  nucleus  is  provided  since  infinite  sets  of  parameters  can  produce  the  same  recovery 
behavior. 
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In  the  recovery  experiments,  when  the  cell  is  aspirated  into  a  pipet  and  held  for  a  short 
time  before  release,  the  nucleus  may  not  have  the  time  to  deform  to  a  significant  extent,  re- 
sulting upon  release  of  the  cell,  in  the  observed  rapid  initial  recoil.  However,  when  the  hold- 
ing time  is  long  enough  for  the  nucleus  to  respond  and  deform  significantly,  a  nearly  Newto- 
nian recovery  is  possible  as  seen  from  Figure  5.6(b).  The  description  of  the  cell  as  a 
homogeneous  non-Newtonian  fluid  ascribes  this  behavior  to  the  fading  elastic  memory  due 
to  relaxation  of  the  Maxwell  element  (Tran-Son-Tay  et  al.,  1991). 

The  reported  finding  that  the  cell  apparent  viscosity  decreases  when  the  cell  recovers 
from  smaller  deformations  (Hochmuth  et  al.,  1993),  or  increases  when  the  cell  undergoes 
larger  deformations  (Evans  and  Yeung,  1989)  may  be  attributed  to  the  participation  of  the 
nucleus  in  the  deformation  process.  In  the  small  deformation  experiments,  the  cell  is  aspi- 
rated into  a  relatively  large  pipet  and  the  nucleus  does  not  have  to  deform  very  much.  As 
the  cell  is  expelled  from  the  pipet,  the  recovery  is  governed  mainly  by  the  cytoplasmic  layer 
producing  a  lower  cell  apparent  viscosity.  In  the  large  deformation  experiments,  the  cell  is 
aspirated  into  a  relatively  small  pipet  and  the  highly  viscous  nucleus  must  deform  leading 
to  a  higher  cell  apparent  viscosity. 

Thus,  the  present  3-layer  model  study  provides  a  consistent  framework  for  describing 
the  variety  of  behaviors,  Newtonian  as  well  as  non-Newtonian,  exhibited  by  the  cells.  Fu- 
ture efforts  will  be  directed  toward  substantiating  this  heuristic  analysis  by  performing  ex- 
periments and  computations  of  leukocyte  dynamics  in  micropipets. 


CHAPTER  6 
SUMMARY  AND  FUTURE  RESEARCH 

6.1.  Conclusion 

In  this  work,  a  three-layer  model  for  leukocytes,  comprised  of  the  cortex,  cytoplasm 
and  nucleus,  is  investigated  by  deforming  the  cell  and  following  its  recovery  characteristics. 
The  hydrodynamics  of  the  deformation  and  recovery  processes  is  analyzed  with  a  mixed  Eul- 
erian-Lagrangian  method  using  the  immersed  boundary  technique.  The  results  for  the  de- 
formation of  compound  drops  under  a  uniaxial  elongational  flow  confirm  those  of  Stone  and 
Leal  (1990).  It  is  found  that  the  inclusion  of  the  core  increase  the  overall  stiffness  of  the  com- 
pound drop,  even  though  its  presence  produces  a  continuous  stretching  mode  of  the  com- 
pound drop  at  smaller  Cacr  than  for  a  single  phase  drop.  Most  importantly,  the  time-history 
of  deformation  is  stored  in  the  compound  drop  in  terms  of  the  shapes  of  its  components. 

The  hydrodynamics  of  the  recovery  process  provides  the  following  insights  into  the  me- 
chanical behavior  of  the  leukocytes: 

(1).  The  cell  apparent  viscosity  can  be  varied  without  recourse  to  non-Newtonian  fluid 
models,  since  the  work  done  in  deforming  the  cell  is  not  entirely  dissipated  by  viscosity 
in  the  cell,  but  can  also  be  stored  in  the  form  of  surface  energy  which  changes  as  the 
shape  of  the  cell  changes.  Thus,  the  resistance  to  deformation  is  not  representative 
merely  of  the  changing  viscosity  of  the  cell,  but  is  also  dependent  on  the  surface  curva- 
tures assumed  by  the  deforming  cell.  This  fact  is  not  accounted  for  in  simple  models 
to  assess  leukocyte  mechanical  properties  (Hochmuth  et  al.,  1993),  leading  to  changes 
in  the  apparent  viscosities  being  attributed  to  the  dependence  of  fluid  viscosity  on  de- 
formation. 
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(2)  .  In  order  for  the  cell  to  behave  as  a  Newtonian  single  drop,  the  nucleus  must  be  suffi- 

ciently deformed  and  the  time  scale  of  recovery  of  the  nucleus  must  be  compatible  with 
that  of  the  cytoplasm-cortex  complex. 

(3)  .  Disparity  in  the  time  scales  between  the  nucleus  and  cell  results  in  a  rapid  recovery  in 

the  initial  stages  which  is  representative  of  the  cytoplasmic  flow,  followed  by  a  slower 
stage  when  the  entire  cell  recovers  as  a  whole.  This  behavior  is  also  observed  for  cells 
with  a  non-deformed  nucleus.  The  non-Newtonian  models  used  hitherto,  such  as  Max- 
well and  power-law  models  do  not  account  for  this  time-scale  variation  due  to  the  inter- 
nal structure.  Thus,  relating  the  apparent  viscosity  of  the  cell  to  the  strain  rate,  as  these 
models  do,  does  not  carry  sufficient  information  on  the  relaxation  times. 

(4)  .  Finally,  assigning  an  apparent  viscosity  to  the  leukocyte  can  be  misleading  without  in- 

formation on  the  shape  of  the  cell  under  deformation  and  recovery  and  can  lead  to  incor- 
rect estimates  of  leukocyte  material  properties. 

In  conclusion,  findings  from  the  present  three-layer  model  qualitatively  explain  the 
physics  behind  the  complex  behavior  of  the  leukocytes  and  reconcile  the  apparently  incon- 
sistent experimental  observations  reported  in  the  literature.  This  work  helps  to  fill  the  gap 
in  the  present  understanding  regarding  the  rheological  behavior  of  leukocytes  and  their  inter- 
nal structural  characteristics. 

6.2.  Future  Research 

In  this  study,  a  computational  methodology  is  applied  to  investigate  the  deformation 
and  recovery  dynamics  of  leukocytes  using  a  multilayer  Newtonian  fluid  model.  However, 
the  potential  of  this  numerical  algorithm  has  not  been  fully  exploited.  The  importance  of 
the  leukocyte's  nucleus  on  the  cell  rheological  behavior  has  not  been  completely  elucidated. 
Additional  work  is  required  to  solve  these  deficiencies.  Some  suggestions  are  given  below: 
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6.2.1.  Experimental  Work 

It  is  imperative  to  develop  a  technique  for  visualizing  and  obtaining  the  necessary  in- 
formation about  the  nucleus  shape  under  different  flow  conditions  in  order  to  provide  a  quan- 
titative description  of  the  mechanical  properties  of  leukocytes.  This  technique,  such  as  the 
fluorescent  staining  procedure,  should  be  incorporated  with  the  micropipet  setup  to  monitor 
the  nucleus  motion  during  the  aspiration  and  recovery  processes. 

6.2.2.  Numerical  Algorithm 

There  are  several  aspects  of  the  code  that  can  be  improved  in  order  to  enhance  the  per- 
formance of  the  present  computational  method.  First  of  all,  the  multigrid  technique  can  be 
implemented  to  accelerate  the  convergence  rate  of  the  code  by  a  substantial  factor  since  the 
convergence  for  the  low  Reynolds  number  flows  usually  are  quite  slow  due  to  diffusion  dom- 
inated process. 

Secondly,  the  lubrication  theory  should  be  added  into  the  code  in  order  to  handle  the 
simulation  of  the  aspiration  process  when  the  gap  between  the  cell  membrane  and  the  pipet 
surface  becomes  extremely  small. 

Third,  if  a  viscoelastic  fluid  representation  is  required  in  order  to  explain  the  dynamics 
of  leukocytes  and  to  be  physiologically  consistent  with  the  properties  of  the  cell  components, 
it  can  be  incorporated  in  a  straightforward  manner  into  the  solution  procedure,  as  done  by 
Hu  &  Joseph  (1990),  and  Yoo  &  Na  (1991).  For  example,  a  Maxwell  fluid  or  viscoelastic 
solid  model  may  be  needed  to  describe  the  cytoplasmic  fluid  or  nucleus  (Dong  et  al.,  1988; 
Dong  et  al.,  1991;  Tran-Son-Tay  et  al.,  1994b). 

Finally,  it  is  essential  to  extend  the  code  to  full  3-dimensional  algorithm  to  consider 
certain  situations  that  cannot  be  treated  under  the  assumption  of  axisymmetry.  For  example, 
in  the  case  of  neutrophils,  where  the  nucleus  of  the  cell  is  segmented,  a  full  3-dimensional 
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description  of  the  cell  and  nucleus  becomes  necessary  in  order  to  obtain  a  faithful  representa- 
tion of  the  cellular  morphology. 
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